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ABSTRACT 


This rsssarch used Huygsns-Frssnst wive optics computer simulations to 
investigate the effects of turbulence strength and inner sarie on the 
normaHzed irradiance variance and coherence length of electromagnetic waves 
pr(H>dg«fing throunh a htfbulent atmosphere. These mvestigations (feveloped 
several fyMeHri,. .or udity of propagation simidations employtoig a numerical, 
split-step, Huygens-r /esr«ei method, end within these gukJelmes. oonsidwed five 
types of turbulence spectrum tnrm scale; (1) zero inrier scale. (2) Gaussian irmer 
scale, (3) Hiirs arKf (4) FrehUch's visoous-convective erihancement inner scales, 
and (5) turbuienoe spectrum truncation from the discrete grid representation. The 
simi^ation results showed that the normalized irradiance variarx^ generally 
decreased <'•30%) below the zero inner scale vaNies in the Rytov regirTte with 
increasir>g inner scale size, but increased monotonicaly in the saturation regime, 
arnf agreed within 2% of the Rytov-T^arsW predictions at tow turbulence strengths. 
The E-flekS coherence length in a spatiaty cortfined beam, with either spherical or 
ptwie wave divergerKe arxf zero irvier scale, f oiowed the Rytov-Tatarski-Fried 
predicitons in die Rytov regime, but deparfed from the theory in the sabiration 
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I. INTRODUCTION 


Turbulwica in a ^nHiflad fluid cautaa random inhomoganaiflas in 
tamparatura and tha indax of rafracflon that scattar and diffract a wava 
propagating throu(^ such a madkim. Analytical parturt>ation tachroquas 
davalopad ovar tha past thraa dacadaa cannot account (Or tha variations 
obsarvad axparfrnantally in tha ampWuda and phasa dMrttxjtiorM of a 
profMgatkig wava whan turtulanoa lavals ara high arxUor propagation paths are 
long Ovar tha past fl fl a a n yaars. numart^ wava optics codas baaad on tha 
Huygar^Frawtal prindpia hava baan dav a lopad to addrass thasa situations 
This rasaarch axiandad thasa codas to includa innar scala in tha spacbum of 
rafTactiva indax fluctuations and to axamina tfta ooharanoa langih and tha 
affacts of innar scala in the Rytov ragima (low turtMianoa stra ngt hs andfor shod 
propagation paths) and in tha saturation ragima (high turtulanoa strengths 
arHl/br long propagation paths) 

As a airtar iWva anatvifai and *^ « * i* *«* davalQoad auktoHnas for 

validity of oomputar simulalions employing Huygsns-f rasnsi propagation ovar 
muflipla slaps (spin-slap method) Thasa g u idsins s includa 
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• For an NxN grid and propagation diatanoa L. the grid elemant size should 

be 

• The maxirnum turbulerioe strerigths C.’and propagation distances L for 
spherical wave propagation with an NxN grid should satisfy 

^ 0.1 N« where P« • 0.497 C.* k"* L"'* . 

• Equivalentiy. the E*fleld coherence ler^th r, should satisfy 

r, 2 2.5 Ax. where r, represents Fried’s coherence length (for 
Sf^MKical waves) 

• Use 2 30 phase screensfsieps for each propagation 

• Use 2 30 propagation realizations to get represerrtattive statistics 

• Phase saeens require low spatial fraquericy correction to gain 5% 
jccuracy in normalized irradianoe variance and as ntuch as 30% in 
coherence length 

• Half Width at half maximum of the irtmospheric MTF erto an iterative fit r, 
provide the most stable p a r am eto rtzatt o ns of the E*field coherence length 

• Telltale signs of aliasing incfude a fine-grained kradtonoe pattern, a boxed 
perimeter to the kradiance patlam. and peaking of the energy toward the 
cantor of the comput^wn grid 

These investigations also examined four choloes of E-ftold sotoce function and 
refined the methods of normaized krad to nc e variance and E-ftold coherence 
length calakation 

These simutottons tocorporatod Ave types of lurbulenoe spectrum kmer 
scale (1) zero inner scale. (2) Gaussian toner scale, p) HVs and (4) FrehliCh's 
viscous-convective enhancement toner scales, and (5) turpuienoe spectrum 
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truncation from the cMsaete grid representation. For the more physicaity 
realistic viscouS'ConvecUve enharwement scale, the computer simulations 
provided the foHovving results; 

• The Hitt «id Frehllch pman^erizations performed almost identicaiiy. 
giving less than 3% difference in normaitzed irradlwioe variance over the 
Rytov and saU^tton regttnes. 

• In the Rytov regime, ttw normalized irradiance variance for an 
approximately spherically diverging beam increased (~30%) and then 
decreased («30%) compared to the zero inner scale values as the inner 
scale size ifKreased. and the E-fleld coherence length r^ decreased 
slightiy (~5%) compared to the zero inner scale coherence length. 

• In the saturation regime, tocreasmg inner scale size gave morK)tonicaiiy 
increasmg rMrmattzed irradiarKe variarKe artol increased the coherence 
length (up to ^50%) compared to ttMi zero inner scale case. 

The effect of the Gaussian toner scale on normalized irradiance variance ar>d 

cohererxto length was also tovestigatod 

These investigations examined the behavior ot the E-field coherertoe 

length for E-fields (beams) thiM were constrained to laleral extent to the grid 

size but whose divergenoe approxtourted that of spherical waves, plane waves, 

and intermediate beam waves The con^M^ stoniattons showed that: 

• For the Rylov regime, toe sph erical and plane wave cMes gave 
coherence len^^ witoto 5% of toe Rytov-Tatarslti-Fried pradictio n s. 

• In the saturation regime, toe E-ltald coherence lengtos for toe spherical 
wave approximalion dea eas e d *25% below toe theory, whtte the 
osherance lertgtos for toe plane wave approxtouMton varied withto 
•5%/^15% of the toeory 
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• The E-field coherence lengths for the beam wave cases were spaced 
between the spherical and plane wave values for the Rytov regime, but 
increased toward the spherical wave values in the saturation regime. 

The spatially confined beams used in these simulations showed a departure in 

behavior of the E-fieid coherence length in the saturation regime from the 

predictions of first order perturbation theory. 

The organization of this dissertation generally follows the preceding 

summary of major points. Chapter II summarizes the theory of wave optics and 

the impiert antation of the Huygens-Fresnei propagation in a computer 

simulation, it then provides an overview of the Rytov-Tatarski theory for the 

effect of atmospheric turbulence, including inner scale, on the normalized 

irradiance variance and E-fieid coherence lengtti as applied to spherical wave 

propagation. Chapter III discusses the detailed choices of physical and 

simulation parameters and develops the guidelines to ensure validity of the 

simulation. Chapter IV presents the computer simulation results, first for the 

inner scale effects on normalized irradiance variance, then for the behavior of 

the E-fiekj coherence length in the Rytov and saturation regimes, and for the 

effect of inner scale on coherefioe length. 
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II. BACKGROUND 


A. CHAPTER OVERVIEW 

This chapter summarizes first the theory of electromagnetic wave 
propagation and the Huygens-Fresnel solution to the scalar wave equation and 
recasts the Huygens-Fresnel solution into a form utilizing fast Fourier 
transforms easily implemented in a computer simulation. It then summarizes 
the Rytov-Tatarski solution to the scalar wave equation which used first order 
perturbation theory and statistical techniques to describe the spatial variations 
of *he irradiance of an electromagnetic wave propagated through a turbulent 
medium. The Kolmogorov spectrum of refractive index fluctuations is 
introduced as well as five types of inner scale which modify the high spatial 
frequency portion of this spectrum. The chapter concludes by describing 
Fried's method for parameterizing the coherence length of the E-field for a 
spherically-diverging wave propagating through a turbulent medium. 

a WAVE PROPAGATION 

Maxwell's equations describe of the propagation of electromagnetic 
waves forough a turbulent medium. Assuming a locally homogeneous, 
isotropic, and linear medium. 




V -ci? • p. 

CD 

V • - 0. 

(2) 


(2) 

d. ,i. LStil, 

9t 

CD 


wfMT* h«tt«d (*} quantttiM mpnatti v^dort Al th« frtqusnciM and ft«M 
•trangths of toitarast. tha madium may ba asaumad to hava zaro iocai fraa 
charga dansity p and zaro (or nagiigibia) oorxSucttvtty o Laaar baama in lha 
atmoapbara aatiafy thaaa oonditioria ExparyJing lha lafl-hand aida of Eq (1) 
arKf comb^)^ fha curl of Eq (3) and tfia dma dartvadva of Eq (4) givaa a 
vactor wava aquabon for tha E-flaW 

v*^ ♦ V(# • V In €) - i» f ^ • 0. (•> 

af* 

Tha mkkfla farm rapraaanta dapoiahzalion of lha E ftatd and ia naglig«>ly amaM 
for propagation through a lurtuidit Mmoaphara (Tatardd, 1961, Lawranoa and 
Strohbahn. 1970) Naglactir>g tha middto M»rm aimpMaa tha vador wava 
aquabon to 

. 0 . (•) 

af» 
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!>)• Fr»tn«(>Kirchofr ctiffraction theory provKfM an approidmM tcily 
solution to this vector wove equMon Foioeting HecM (1987). the ip heric e l 
wi^ solution m 

Eiu ). •-« ->. m 

Where E, is a constant and k ■ 2s/X Using Eq (7) w(0) lha time dependence 
factored out. the Kirchoff integral theorem becomes 

Which must be evalueted over a surface S endosing the field poM P. Figure 
(1) Wiartfatei the relettonship bet wee n the dtstarKses ^ and n if the w a v e l e r H Uh 
X « (. n . Eq (8) reduces to Pie Fresnet-Kiithofr dWir acU on ftirmula 

• T-' //.^ '"*> ** *** 

Where K(6) represents the obliquity factor 

i nte g r a tin g Eq (9) over the half sphere S (shown in Fig. (1)) as the 
radius approaches infinity, the electric tieid amplitude » appreciable only over a 
Ande aperii^ area A of the hemisphere's 1UM side Then Eq (9) reduces to the 
more famitiar Huygeris^Fresriel expression 
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PIpM* 1 F rwn t t Wrc ho ff dWfr acO on g totn t ry showing rsiottonship b tws n 
souros. sporturs A. surfisoo 8. and fItW point P. 

A ^ 

Ths ootnpm aparlurt hmdion oonlains tha aouroa aphahcal wava factor, 
a]ip(itOi^ T?w E-flaM at a point P it a function of tha E-fia(d ovar a finila 
apartura A. 
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The Huygene>Fresnei formulation . Eq. (10), serves as foe basis for 
calculating foe propagation of a wave through turbulenoe. This formula can 
now be cast foto a form easily implemented in a propagation computer 
simulation. Applying foe Cartesian coordinate system of Fig. (2), foe Huygens- 
Fresnel prfodpie becomes. 

eVA • T //. -- j-■ 'fW 

^ ^i* *\t- pi* 



Aperture 


Flgurs 2 Propagation coordfoate system showing relations between foe 
aperture variabfos and foe field variables. 
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^ reprtstnis the position vector of the observation point P in the x-y plane and 
p represents the position vector of the radiating point in the aperture plane. 

The paraxial approximation assumes that 

\f\*z 

where z measures the lortgitiKflnai distar>oa of the point P from the aperttfe. 
Consequently, the obliquity factor K(OH1'*^oa0)/2 becomAs approximately 1 
and the denominator of the integrarKt is approximately z .x>dman.1968), 

EVA dA. ii») 


Fottowing Roberts (1986), the Fresnel approximation assumes that the lateral 
displacement between the radiating point and the observation point 
I p - p I « z. Then, using the lowest order terms for the series expwision 
of tfM complex expor)entii^s argument, Eq. (13) becomes 

. - -itefji - |iZf 

EVA • E{pSt) • ‘ * ' 


XI 






(14) 


. a«y 

' f; ‘ I* 
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The fBtilor exp<-2xiz/X.) represents the change in phase of the optical wave from 
the center of the aperture plane to the plane of the observation point P. It 
applies to the v^e E-fieid and deperuls only on the propagation distance z. It 
does not affect the phase variation across the E-field nor the amplitude of the 
E-fiekj. and will be dropped in the following expressions for simplicity. 

Equation (14) can be recast into a form utittzing Fourier transforms. 
Expandirni the aperture field E(^0) in Eq. (14) using the Fourier transform 
identity 

£(p.O) . E(p',0), (1«) 

Eq. (14) becomes 

E(f,2)- ~ fdp' E(p',0) ] <“*•) 

(Note; these equations use spatial frequency / (m ') instead of spatial 

wavenumber k > 2 s / (rad/m) because the discrete Fourier transform 
computer subroutines used in these investigations are written with spatial 
frequency f .) Making the change of variables f'mf-p gives 
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E{t^ - J(-dr) [ /dip' ] • 


—1^1* 


( 17 ) 


fdf9*'^fdp' • E(p'.0) [ fdt' 9 


-m*y ^ XX 


i* iVi • 




]■ 


The last intagiation over it the Fourier transform of a Gaussian function, 


/ 


df/ a*' 


-1^1* 


/Xr 


( 18 ) 


Substituting pfor p', the E-fiekf of Eq. (17) becomes 

E(A^) « Jdf a '-Mfl* /up £(p,o) #-*■'>. (1*) 

Equation (19) may be symbolicaNy written as 

E(A2) • IFTl EriE(p,0)]I, (2®) 

where FT and IFT represent the two-d^ensionai Fourier transform and inverse 
Fourier transform, respectively. Equation (20) expresses ttie E-fiefd at a 
propagation cHstanoe z in terms of Fourier transforms of the E>fleld at z » 0 
represented in a Cartesian coordinate system. This form of the Huygens- 
Fresnel principle is useful for a emulation that subdivides the propagation path 
into a sequence of short Fresnel propagation steps. 
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C. EFFECTS OF TURBULENCE ON THE PROPAGATED WAVE 

1. Imdhinot VarltiiM and Innar Scala 

The next atop in modeling the propagation of an E-field through a 
turbulent medium requirea a atatiatical characterization of the turbulence and its 
effect on the apatiai atobatica of the propagating wave. To begin understanding 
the effect of turbulence upon the wave, consider a spherical wave incident upon 
a merbum rarxiomly varyirni index of refraction, as shown in Fig. (3). For 
amaH variations in the index of refraction about the mean value, scattering 
occurs predominantiy in the forward direction (TatarskI, 1961). The regions of 
index of refraction fluctuation accelerato/retard portions of the wavefront over 
short propagation diatonoes and cause variations in phase across the reference 
spherical wavefront. Subsequent (flffraction and interference then create 
variations in irradianca across the wavefront. Compared to propagation through 
zero turbulence (Fig. (4)), a sphericaRy diverging beam now exhibits significant 
variations in phase stkI amplitude (irradianca) (Fig. (5)). 

A statisticai description of the refractive index fluctuations is 
needed to calculate the E4ield variations. The index of refraction, n, in air 
depends upon the temperature T and pressure P (Tatorski, 1961). At visible 
wavelengths, 
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Flgurt 3 Wav* propagating through a madium having random iNK>mog6naitia8 
in indax of rafraction. 


n - 1 - 72x10^ (21) 

T{Ki 

Tha tamparatura proflia in tha atmosphara poaaassas signincant atratification. 
as shown in Fig. (6) (Waitars, 1994). Valocity shaar batwaan difTarant layers in 
tha atmosphara causes turbulanoa which disrupts tha intarfaoa batwaan these 
layers, mixing regions of dNTarant tamparaturas. Tha valocity fluctuations 
ganaralad by tha turbulanoa at tha intartaoa batwaan layers ganaraily foNow tha 
Kolmogorov spectrum. 
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Rgura 4 Irradiance plot for a spharfoa^ diverging beam propagated through 
zero turbulence. 






ngurof Irradiano* ptol for sph«rfcaly diverging bMm propagate 
turtouiafKa. 







Ftgui* • Vertical proWa of Itmparaliirt in the atonotpharo. 

♦<«) - 

aaauming iaolropic and homoganaoua turtuianoa k rapraaanta tpattal 
wavanum{>ar in rad^ Tamparalura iaa Ipaaaiva additiva* (Tatardd. 1961) in 
tha atinoapfMca. maaning that K doaa not affaol tha dynantica of lha turbulanoa. 
Thua. tha thraa-dimanaional tpadrum of tawparatura fluctuationa alao fotioart 
tha Kolmogorov apactnjm SInca Eq. (21) raMaa tha indax of fafr a ctio n to 
lamparaiura, «ia ipacarum or ranracava moax aucaiaDona la aiao isovnogorov 
(TatanW. 1661) 

♦.(x,/) - oon qfti) x*”^. (») 
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C,’(z) is a cx>nsUint of propwtionMMy indicating the «drengtti of the 
tuftHiience. in reality, this spectrum only appUes to an intermediate range of 
wavenumbers known as die inertia subrange, shown in Fig. (7). The 
boundaries at the high and low spatial wavenumbers are the inner and outer 
scales, respectively, and are cNscussed below. 



Figure 7 Atmospheric spectrum of retractive index fluctuations showing inertiiy 
subrange and outer and Inner scales 

this u n derst a nding of the effect of t u rbulence on a 
propagating E-fleld. TatarsM (1961) derived expressions for the statistical 
p r o p arbas of the ampWude. phase, and irr a d la no e of a spherical wave 
propagatlrig through a turbulent medkan as foiowt Consider a scalar 
c om ponent of the vector weve equation Eq (6) and assume a hamtonic time 
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I 

dependence for the E-field 

Em - uin 

Perform the time derivatives and introduce the wavenumber k = 2iUX and the 
irKtex of refraction 

/i = Cv^ <26) 

(where c » velocity of light in vacuum, p » magnetic permeability of the 

medium, e - electric permittivity of the medium ) to get the scalar wave 
equation 

0. (26) 

Expressing the index of refraction as 

n{t) * 1 » fiXt), (27) 

where ln,(^) I « 1, and usirtg the Rytov metfrod of smooth perturbations that 
employs u * exp(H0. 4' * V, T, + ..., and u = u^ ♦ u, + ..., the wave 
equation becomes 

♦ 2ilr*n,(/) - 0. (28) 

This has a sokrtion 

▼i(') - fhin dV'. (29) 
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■ 




for a tpharical wav* 

u^\f\) « (O . oontfanfi (30) 

In 


and for small wavaiengtht X € whara ^ (tha innar scaia) charactarizas tha 
siza of tha smaHast fkictuatiorw in tha tndax of rafraction. tha FrawMN 
approximation in Cartasian ooordinatas raduoas Eq. (29) to 




aKp( Mr 


zHx^* ♦ (**♦/*) - 2zz* 

__ 

z'(z-z^ 



Tha statistics of tha turtHjIanca appaar in H'f and n, in tarms of spactral 
axpanskxts whara n, oorUairts tha thraa*dimansional spatial fraquancy spactrum 
of rafractiva irtdax fluctuations, 0 ,(k). Tha Ryfov approximation introduoas tha 
log ampWuda fhjctuafion X 

X . h 4. (>») 


whara A and A, are tha ampiitudas of tha turt)ular>ca paftuft>ad E-flald u and 
tha firaa spaca E*fiald u. 


u • Ao** 

■ Ap 

Than, after an axtendad sarias of man^HJlations (saa Tatarski, 1961) that 


(M) 


aasuma loctf homoganaity arxj isotropy, tha varianca of tha log amplituda 
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fluctuatkm X of a spharical wave propagating over a distance L becomes 


? . 4**** IJdK > ^ ♦.M) 

Where k • 2tUX, and k « spatial wavenumber (rad/m), which is used here 
instead of spatM ftequericy f (1/m) to be consistent with previous work. 
Assuming that the irradianoe follows a log normal distribution (Tatarski, 1961), 
the normalized irradiance variance is a function of the log amplitude variance 

- «p( 4 ? ) - 1. 

The assumption of isotropic. homoger>eous turbuierKe gives a 
Kolmogorov spectrum of refractive irfoex fluctuations of the form 
^.(k.z) s 0.033 C.’(z) k'"'’ (Tatarski. 1961). Assuming, for computational 
convenience, that this k'"^ power law dependence holds over the entire range 
of spatial wavenumbers k. and that the turbulence strength C.’(z) is uniform 
along the path, (hen integration of Eq. (34) gives the log amplitude variance for 
a spherical wave as 

X* - 0.124 Cj (“» 

The normalized irradianoe variance of Eq. (35) becomes 
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of 


7 * 


«(p(4 ?) - 1 « mpC 0^ Cj - 1. 


(37) 


Low turbutonce (smaN and/or short path lengths make the exponent small 
enough to apply the approximation exp(a)« 1*Kx. giving 

^ - 0.407 Cf - pI , (33) 

where the new parameter pj, defirted by this equation, serves as the baseline 
normalized irradUmoe variance for comparing modifications to the spectrum of 
refractive if>dex fluctuations, pj also ^iitates plotting normalized irradiance 
variance over a broad range of integrated turbulence strengths. 

The above steps characterize the effects of the turbulent medium 
upon the irradiar>oe statistics of a spherical wave assuming a simple 
Koknogorov spectrum of refractive ir>dex fiuctuatiorvi. Incorporating a high 
spatial frequency roHoff at the inr>er scale, as shown in Fig. (7). modifies ttre 
kradiance statistics. The inner scale ^ represents the physical size where 
viscosity of the medium smooths the velocity fluctuations by dissipating kinetic 
errergy and tttermal diffusion smooths the temperature fluctuations (thus 
refractive index fiuctuattons), removing the turbulent character of the medium at 
this scale. This smoothing causes a roitoff of the high spatial frequency energy 
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of the refractive index spectnim. Flatty, Wang, and Martin (1993) write the 
three^imensional isotropic spectrum as 


•(k,^) » 0.033 F{kI^), (5») 

where F(k^ represents a particular functional form of the inner scale. The 
parameter consists of the spatial wavenumber k (radians / meter) and the 
inner scale parameter (meters). For zero inner scale (i.e. no high spatial 
frequency rollofr), 

= 1. 0 < K < -. (40) 

For theoretical and computational convenience, a Gaussian rolloff inner scale is 
often used (Tatarski. 1961) to represent the high spatial frequency rolloff from 
viscosity 


The more realistic viscous-convective enhancement inner scale advocated by 
Hill (1978) and Frehlich (1992) exhibits an enhanced spectrum for the 
waverujmbers slightly less than the rolloff point. Viscosity attenuates the kinetic 
energy and lowers the velocity fluctuations over regions near the inner scale 
size and smaller before thermal diffusion can smooth all the temperature 
fluctuations within tiiese regions. This disparity alters the temperature and 
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i 

corrasponding index of refraction specfra. Inner-8cale>8ized patches of air have 
different temperatures but have little internal velocity variation. Since the 
Kolmogorov spectrum is based upon the spectrum of velocity fluctuations, these 
residual temperature fluctuations cause an enhancement to the Kolmogorov 
spectrum around the inner scale size. At higher spatial frequencies beyond the 
inner scale, thermal diffusion does smooth out these residual temperature 
fluctuations and the spectrum falls off sharply. 

Hill provides a plot of F(k^ versus for the viscous-convective 
enhancement, that is suitable for implementing the viscous-convective 
enhancement inner scale for numerical integrations and computer simulations 
(Flattd, Wang, and Martin, 1993). Frehiich presents a similar characterization of 
the viscous-convective enhancement inner scale based upon laser scintillation 
measurements and provides a four parameter fit to describe this version of the 
viscous-convective enhancement inner scale. 

Since most numerical simulations, such as the one used here, 
utilize discrete Fourier transforms, the spatial frequency grid mesh chosen for 
calculations introduces a maximum ^atial frequency established by the 
Nyquist sampling criterion (discussed in Chapter III). This limit creates a grid 
cutoff inner scale at k,^. 




0 < K < K, 
K > K, 


(42) 
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Figure (8) illustrates these five inner scales for eight inner scale 
sizes. Solid curves represent the Hilt ver»on of the viscous-convective inner 
scale, dashed curves represent the Frehlich version of the viscous-convective 
inner scale, and dotted curves represent Gaussian inner scales. The solid line 
with the box shape represents the numerical grid cutoff inner scale at = 
318 rad/m for a 1024x1024 mesh. The inner scaie values of 2, 3, 4, 5, 6, 7, 
10, and 15 cm refer to stratospheric propagation over 200 km, assuming an 
optical wavelength of 500 nm. 

The outer scale l^, corresponding to the spatiai wavenumber 
- 2ii/Lo, represents the upper size limit to which the Kolmogorov spectrum 
applies. For scale sizes larger than Lq (or spatial wavenumbers less than k,„), 
the refractive Index fluctuations level off to a finite value as k approaches zero. 
Physically, such a limitation exists as the spatial wavenumbers approach zero 
because the refractive index fluctuations cannot become arbitrarily large, or 
equivalently, the energy represented by the spectrum must remain finite 
(Tatarski, 1961). Estimates of the outer scale for the stratosphere lie in the 
range of tens to hundreds of meters. However, the outer scale proves to be 
much more problematic to include in the spectrum of refractive index 
fluctuations because the atmosphere is anisotropic at these large sizes 
(Tatarski, 1961). The Von Kdrmdn spectrum (Tatarski, 1961), 
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Wavenumber, K (rod/m) 


FIgurt 8 Inner scale function F(k^ the grid cutoff (box). Gaussian 
(ck>tted). Hill (solid), and Frehlich (dashed) inner scales, plotted for 1024x1024 
grid, L s 200 km. and X » 500 nm. 






•» 


0.033 Cj 

* I*-”'* * 


( 43 ) 


which has been used to Incorpoiate an outer scale for some calculations, is 
based upon computational convenience more than physical understanding. 
Additionally the computer simulation results were compared against the Rytov- 
Tatarski predictions, which were derived without an outer scale. Therefore, 
these investigations did not incorporate an outer scale. 

2. Atmoapherfc MTF and Coheience Length 

Propagation through a turbulent medium not only affects the 
irradiance statistics, but also reduces the spatial coherence of the wave as 
characterized by the transverse coherence length. Fried (1966, p. 1372-1379) 
derived a long exposure modulation transfer function (MTF) for a spherical 
wave propagating through a turbulent atmosphere and used this to 
parameterize the E-field transverse coherence length, r^. Following the method 
and notation of Fried (which used spatial frequency f in ( m'^)). consider the 
spatial Fourier transform of the intensity of an E-field propagated through the 
turbulent medium and imaged by a thin lens 

x(f) • Bf at U(J0 (44) 

where u(j^) represents the E-fleld in foe image plane and B represents a 


27 



normalization constant. Fried caiis this the "MTF of the image-forming optical 
system**, assuming a unit impulse (point) illumination and anticipating the fact 

that the ensemble average of t(/). which incorporates the effects of 

atmospheric turbulence, turns out to be real. Utilizing the Fourier transform 
property of a thin lens, 

uiJi) ^Afd9 Ui^) 

where A is another normalization constant and U(^) represents the E-fieid in 
the plane of the thin lens aperture, the MTF becomes 

t(/) ^ A^Bf d» - Xf^) U{^). (^) 

Expressing the E-field U( P) as the product of a zero turbulence propagation 
part, W(^), and an atmosphere-induced perturbation, V(^), 

U(9) - W{9) V^(P) - W{^ 

W(represents the uniformly illuminated aperture function. I{9) represents 
fluctuations in the log amplitude of the E-field and ^ represents phase 
fluctuations. 

Taking the ensemble average <> of Eq. (46) over many 
realizations and substituting in Eq. (47) gives the iong exposure MTF 
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<^ih)iongmpo$um = V^(P) >. W 

Fried denotes the expectation of the fluctuations as the "atmosphere's MTP, 

MTF^f) = < V{^ - Xflf) V{9) \ ( 4 «) 

and Hufhagel and Stanley (1964) call this the average mutual coherence factor. 
Substituting from Eq. (47). 

( ^*(9 - XRf) V{9) ) = \ (#0) 

Using the fact that ^ and I are Gaussian random variables. Fried shows that 
the atmosphere's MTF reduces to 

• . -4 (61) 

< - XRf) V{V)) ^ 9 ^ ' 


where D(XR |/|) represents the wave structure function and is related to the 
phase structure function and the log amplitude structure function D, by 

D(/) » D^{l) ♦ when r* 19- XRf\. («) 


The log amplitude structure function D, and phase structure function are 
defined as 


o,W ■ (|l(^ - ), 

o«W ■ (l*W - ♦(»-x/V)i* >. 


(SS) 
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•nd cteicribe th« expactad variation in tha log ampHtuda and phasa. 
ratpactivaly. at Uto points a distanca XRf apart. Danoting tfia MTF of tha 
apartura function W(as 

t,(/) ■ A*B f W(P), (M) 


and using tha tact that tha atmospharic MTF of Eq. (51) is indapandant of P, 


tha long axposura MTF of Eq. (48) now bacomas 




ocxaifD 

$ 


(56) 


or. 


-4 o(W 

a * 



( 65 ) 


Tha long axposura MTF is raiatad to tha mutual ooharanca of tha 
E-fiald. Tha mutual ooharanca function (MCF) dascribas tha autocorrelation 
batwaan tha E-fiald at two points and is defined as (Goodman, 1985) 

MCF - ( (/•(/►.I,) >, («T) 

whara U raprasants tha oomplax E-fiald. For these investigations t^^t, and 
explicit time notation will be dropped. Asstxning homogeneity, tha MCF 
becomes tha spatial autocorralation of tha E-fiald 
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UCF{t) - / dt' U*{F) U{r * f). 


(M) 


Taking tha ansetnbie avarage of Eq. (46) shows that tha long axposura MTF 
comas from tha avataga mutual coherence furKtion (autocorrelation) of the 
E-fiaid in tha aperture of tha Ians 

(t(/)> ^ A*B {f - X«f) U(^ ) • (( (/•(©) U{^*XRf))). («•> 


Rewriting Eq. (56), 





< W(t') W{F*/) ) 


m 


Equation (60) is suitcbla fOr investigating tha wave structure function D(r) sir>oa 
computer simulations provida tha E*fialds U and W on tha right hand side of 
Eq. (60). 

Fried (1966, P.1380'1384) derived tha wave structure furKtion for 
a spherical wave based on tha thraa'dimansional spatial frequency spectrum of 
refractiva index fluctuations, ^.(k.z). 

0(1) . 6«**» dt // [1 - ♦.M . dc. (*1) 

where z represents the position along the optical path length, 0 < z < L . Fried 
analytically solves Eq. (61) for a spherical wave and simple Koimogorov 
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furbultnot that hM an indttc of rsiricllon itruclur* ftjnclion 


D,M - cj(z) 1^, (“I 

and a thraa-cHmansioniri spattai tpadnim of rafractiva vviax fluctuationa 

#(k.^ - 0X83 CjW k (•») 

SuhstHut^ Eq. (63) into Eq. (61) and intagrattig ovar spatial fraquancy k givas 
D(/) -ZSi k* dz Ci(z) (^. (64) 


Assummg C.*(z) s constant along tha opticai path givas 

D(/) . 1.060 k* d Lr^. (•») 


This aquation ralatas tha wava stnjctura function 0(r) to physicai paramatars 
k«2ii/X, C.^ and L for tha propagation through turtMianoa. FoAowing Friad arKl 
daflnk)g tha constant r. 


4 ./_ m _ 

'1X80 Ir* Cf i. 



( 66 ) 


Eq. (66) tiaoofnas 

D{/) . 6J6 (-JT- (6^ 

Tha parama^ r, c har actart ia a tha o o haranoa lan(^ of tha E-I)ald bacauM tha 
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resolution ailowod by the turbulent atmosphere is 0/^ • X / Tq (Fried, 1966, 
p. 1380-1384). 

Substituting Eq. (67) into Eq. (60) gives 

((wn uv'*D > > ( 68 ) 

< ivto W{f*r) > ’ 

and now dire^ relates a measure of the coherence length, rg, to the E-fields 
produced by computer simulation. Similarly, substitLting the integral form of the 
wave ^ructure function, Eq. (61). into Eq. (60) and numerically integrating for 
spectra ^.(k,z) with different inner scales allows comparing theory and 
computer simulation. This comparison tadlitates validation of the computer 
simulations incorporating an Inner scale at low turbulence strengths where the 
perturbation<based theory remains valid, and provides a mechanism to explore 
conditions of high turbulence strer>gth and/or long propagation path lengths (the 
saturation regime) where the perturbation theory may no longer hold. 
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III. COMPUTER SIMULATION 


A. PROPAGATION CODE 

These investigations utilized a wave optics computer simulation code 
written by Brent Ellerbroek at the United States Air Force Phillips Laboratory in 
Albuquerque. New Mexico that incorporates phase screen generation routines 
written by Greg Cochrane, also with the Riillips Laboratory. The code, known 
as YAPS (Yet Another Propagation Simulation), is a general purpose adaptive 
optics simulation code written in FORTRAN that models optical propagation 
through a turbulent atmosphere, sensing of the wavefront with Hartmann 
elements (Hudgin, 1977} and a CCO array, optimized phase correction 
calculation, and phase compensation via a deformable mirror, all within a time 
indexed framework. These investigations utilized only the propagation portions 
of the code, modified the code to parameterize turbulence strength with 
instead of rg, added different source configurations, incorporated Gaussian, Hill, 
and Frehlich inner scales, and reduced the memory requirements. 

The computer simulation utilizes the sput-step method to simulate 
propagation of the E-*ield through tijrbulence, as illustrated in Fig. (9). A source 
E-field (left) was propagated in steps (with zero turbulence) over the 
propagation distance L with a random i^se screen applied to the field at each 
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Z - 0 Z “ L 


Figm 9 Computer simulation using the split-step method. The source (left) is 
propagated in steps out to z » L, with phase screens applied at each step. 

step. This method was based upon the extended Huygens-Fresnei principle 
(Yura. 1992) 

£,-jfE a - 5 ^ Af(9) •’ . <•») 

which is Eq. (10) with an extra e'^ factor that incorporates the random variations 
in log amplitude f and phase ^ 

*» . •' •'♦. (70) 

For paraxial propagation, if the distances r\ are short enough, diffraction and 
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inMrfBTonos not havo a chanoa to cauaa significant log ampiitijcle 
fiuctuations f, allowing the approximation 

e» * <7l) 

Inserting this approximation into Eqs. (11) • (19) gives 

E{r,z) = Jdp E(^,0) »/♦« (72) 

Thus, a single step of the split-step mefiiod requires applying the phase screen 
e'* to the E-field and propagating the result a distance z using the Huygens- 
Fresnel propagator. 

To implement this sequence of steps, the YAPS code followed a user- 
defined list of tasks (stored as an input file) and called the appropriate routines 
to accomplish those tasks. For these investigations, a typical list (see the 
Appendix) first set up the initial parameters, including random number seed, 
wavelength, number and size of fields, source characteristics, and phase 
screen size. The list then initialized the field grid and applied the initial source 
distribution, propagated the field in steps, generated and applied phase screens 
between steps, and finally saved the complex values of the propagated field. 

The YAPS propagations were run on Sun SparcStation-IO’s having 128 
megabytes of RAM. This memory consb^int allowed a maximum grid size of 
1024x1024 (choosing N as an integer power of 2) for the current version of 
YAPS. Each run consisted of 30 propagations, each with 32 phase 




screens/steps, and required ~16 hours to complete. These investigations 
normally utilized nine SparcStation-10's simultaneously doing runs with different 
parameters. In all, these investigations consumed roughly 6000 hours of 
computer time. A small Cray mainframe was available with more RAM, but was 
not utilized extensively because the multiple SparcStations provided more 
overall computational capacity. 

Once the YAPS code had generated and saved the realizations of the 
E-field propagated through turbulence, separate routines written in interactive 
Data Language (IDL) analyzed and displayed the fields. Analysis included 
displaying two- and three-dimensional plots of the intensity field, calculating the 
dependence of intensity and normalized irradiance variance on radial distance, 
calculating the normalized irradiance variance over the central portion of the 
field, calculating the atmospheric MTF and corresponding coherence length, 
and calculating Strehl and intensity ratios (defined in Section E below). 

The wave optics computer simulation required many choices to model 
the stratospheric propagation scenario and to ensure validity of the simulation. 
These choices included wavelength, propagation distance, inner scale, grid size 
N, the physical size of each grid element Ax, the maximum strength of 
turbulence the source function of tiie E-field, the number of phase 
screens, low spatial frequency corrections, number of realizations, and methods 
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of calculating average statistical values. The foliowing sections address ttiese 
and other choices and develop guidelines for validity of the simulations. 

a PHYSICAL PARAMETERS 

Wave optics computer simulation of the E-field requires selection of the 
parameters that describe the physical properties of the propagation. Usually a 
specific propagation scenario has been selected for modeling, giving the 
desired wavelength X, wavenumber k=2ic4X, propagation distance L. and the 
range of turbulence strengths C„^. The stratospheric propagation scenario 
chosen here used X = 500 nm, L = 200 km, and = [IxlO"*’, 1x10'*l m"*®. 
From these physical parameters, other useful scaling parameters occur, such 
as the Fresnel wavenumber (Martin and Flatt6, 1988) 

• (i)’" = ^ 



the Rytov-Tatarski normalized irradiance variance (Flatt6, Wang, and Martin, 
1993) 

^ = a Cj . (M) 


where a-1.23 for plane waves and a-0.497 for spherical waves, and Fried's 
coherence length for spherical wave propagation (Fried. 1966, p.1380-1384) 
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(78) 


^0 = 


6.88 


m 


2.91 Ir* dz Ciiz) (!)•' 

= (-(Ay aifisttwf Cj). 

M.09 k^d U 


The inner scale of the turbulence, often is not known exactly, but values can 
be estimated. For stratospheric propagation, estimates of the inner scale are 
around 1 -15 cm (Beland, 1993). Once the inner scale is known, the Hill and 
Frehlich versions of the viscous-convective enhancement inner scales (Hill and 
Clifford, 1978) (Frehlich, 1992) were implemented with the parameter , 
where k represents spatial frequency. The outer scale Lg for stratospheric 
propagation lies in the range of tens to hundreds of meters but, again, due to 
the difficulty in parameterization, was not included in these simulations. 


C. COMPUTATION GRID 

The E-field was represented as an NxN array of complex numbers in the 
plane perpendicular to the axis of propagation . Generally, N was chosen as 
large as possible, consistent with the amount of computer memory available 
and the time required to run a simulation. Larger grid size N allowed a wider 
spatial extent of the field and/or sampling of higher spatial frequencies (denser 
mesh). The physical size that each grid element represented had to be chosen 
to ensure validity of the computer simulation. 
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Knepp (1983) discussed the relation of grid element size Ax and physical 
grid extent NAx to inner and outer scales, to proper sampling of the phase 
screens, to angular spreading of the field, and to proper sampling of the spatial 
frequency quadratic phase factor in the Fourier transform formulation of the 
Huygens-Fresnel propagation. For ^e latter, the quadratic phase factor in 
Eq. (20) takes the form 


♦ 


A a. * 


(78) 


where k is a spatial wavenumber in rad/m. Applying the Nyquist criterion, 
which requires that this quadratic phase change by less than it across one grid 
element Ax, Knepp derives 

^ ^ 2 N i ^ xf ( 77 , 


Roberts (1986) applied similar techniques and derived Eq. (77) without 
the factor of 2. Again applying the Nyquist criterion to the quadratic phase 
factor in the Fourier transform propagator Eq. (20) but utilizing more succinct 
differential notation. 
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and, 


■ AA2 X X L 4hw < 


m 


where represents the spatial frequency with ttte maximum rate of change of 
phase. At the edge of the grid, f,^ is 


2 2 N^x 2 Ax 


(80) 


After substituting and rearranging, 

L < (81) 


When the propagation distance L is known, as for a specific propagation 
scenario, and when the grid dimension, N, is known, then the grid element size 
is 

AX>^. (M) 


Sampling of the phase factor in the spatial frequency domain must meet 
the Nyquist criterion as a minimum. Some circumstances may warrant applying 
even stricter sampling criteria such as restricting the phase to change by less 
than an across one grid element, 0 < a < 1. Furthermore, suppose the 
propagation involves spatial frequencies out to where 0 < p < 1. Then 
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64 < ax. and tha maximum tpatiai firaquancy ia Substituting thasa 
constraints into tha analysis givas 


6x> 


u 

N 


( 63 ) 


Whichavar is usad, Eq. (82) and (63) provide simple formulas for choosing the 
minimum grid alamant size for tha propagation. 

Spherical wave propagation places an additional constraint on the choice 
of grid element size. In the parabolic approximation, a spherical wave has a 
quadratic phase curvature which represents divergence from a point source a 
distance S (focal distance) away 

♦ * 


where p measures the radial distance from the propagation axis. Roberts (1986) 
analyzed the sampling criteria for fois phase just as for the spatial frequency 
phase (though applied to the case of a two-step Fourier transform Huygens- 
Fresnel propagator and not applied specifically to spherical waves). Proceeding 
as above arfo applying the Nyquist criterion to this phase factor 


it. A. (. it. 

tfp dp X S Ap 


( 88 ) 
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( 86 ) 




Ap 2 n < 
A 5 


n. 


represents the maximum radial distance. Using corresponding 
to the edge of the grid and Ap ■ Ax. 


A 

PmtK = -2 = 


Af AX 
2 ‘ 


( 87 ) 


Substituting Eq. (87) into Eq. (86) and rearranging, 

(•») 


Again generalize the analysis by requiring the maximum phase change 
across one element to be an (0 < a < 1), and requiring the field energy to be 
confined within a region of radius yp.^ , (0 < y < 1). This gives 


Ax< 


N 


XT7 

N Y* 


( 89 ) 


Combining Eq. (83) and (89), 

TTe 

N Y* 


HI 

N a 


< Ax< 


( 90 ) 


Spatial frequency sampling considerations for the Huygens-Fresnel propagator 
have placed a lower limit on Ax, while spatial sampling considerations of the 
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quadratic approximation for spharical wava phasa hava placad an uppar limit on 
Ax. 


Tha computar simulations usad in thasa invastigations axaminad high 
lavals of turbulanca that introducad anargy into tha highest spatial frequencies 
reprasantabla on tha grid and that scattered anargy to tha edges of the spatial 
grid. Tha simulations also assumed that tha Nyquist sampling criterion was 
sufficient. These considerations spadfiad tha parameters a = p = Y - and 
lead to 



< Ax< 



( 91 ) 


Additionally, these simulations propagated a spherical wave from the point 
source at the origin (focus) out to a distance S, i.e. S = L, making the 
inequalities in Eq. (91) become the equality 



( 92 ) 


These choices unambiguously determine the guideline for grid element size for 
spherical wave propagation based upon sampling considerations in the spatial 
and spatial frequency domains. Equation (92) also determines the grid element 
size for plane wave propagation since minimizing the grid element size of 
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Eq. (82) optimizes the sampling of high spatial frequency distortions caused by 
turbulence. 

Some simplifications were made in the above analysis. First, the 
maximum spatial distance and maximum spatial frequency used in the 
quadratic phase factors were chosen for the nearest edge of the grid and 
correspond to the radius of the largest drde that will fit inside the grid. These 
choices disregarded the comers of the grid, but this omission should not have 
affected the simulations greatly because the majority of the energy in both the 
spatial and spatial frequency domains was confined within the radius of the 
drde to minimize aliasing. Second, the Nyquist criterion was applied to the 
phase change across the x or y dimension of the element. Analyzing the phase 

change between opposite comers of a grid element increases Af or Ap by v/2 . 

Finally, an NxN grid with N even requires pladng the (x,y) » (0,0) point at a grid 
point, such as (N/2 ■•'1, N/2 ■•■I) that is not the exact geometric center of the 
grid. The distance to the nearest side is (N-2)/2 elements instead of N/2, which 
is a minor change for large N. Incorporating tfiese three considerations into the 
analysis for a spherical wave gives 
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How»v«r, squaring Eq. (93) «k 1 raammging laads to 




i<i_«Us 


♦ (/W-2)‘ 

SO that a spherical wave propagatkKi with these additionai constraints should 


m 


only be carried out over a maximum distance less ttian S/4. Since these 
investigations needed to propagate a spherical wave over the full focus 
distance S starting at the source, the simpler expression Eq. (92) was used to 
determine the grid element size fOr these investigations. As a comparison, 
FlattA, Wang, arni Martin (1903) used a grid eiement size of 0.7mm fOr a 
1024x1024 grid and 0.5mm for a 2048x2048 grid with XL = 0.000638 m^. 


Equation (92) with these parameters prescribes grid element sizes of 0.78 mm 
and 0.55 mm, respectively, which are ~10% larger to meet sampling 
considerations for the Huygens-Fresnel propagator. 

The spW-step method in the ccxnputer simulation divides the optical path 
into multiple steps. However, the dMivice L used to determine the grid 
eiement size must correspond to the total propagation path length and not the 
step size. To justify tNs, consider the vacuum propagation of a field across die 
distance L » Az. If prc^xagation occurs in a sin(^ step, then the Fourier 
transform Huygens-Fresnei propagator, Eq. (20). becomes 


E(Az) - IFT{ a '**^'* F71E(P)]]. 
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If ttw path has two steps, then 


e(Ar) - FTlElizl2)i\, 


m 


and 


E{Azf2) = iFT[ 9 » FT\EiO)]]. 


(97) 


Substituting Eq. (97) into Eq. (96), 


E(Az)-/F71a * FT\/FT[9 * Fr(E(0)llll. ' * 


But ihe middle Fourier transform/ inverse Fourier transform pair cancel each 
other, ar>d this two step vacuum propagation becomes the one step propagation 
of Eq. (95). Thus, the grid element size must be chosen to correspond to the 
total distance Az - L since multi-step vacuum propagations mathematically 
reduce to a single step of L. 

D. SOURCES 

The choice of a given propagation scenario significantly affects the 
irradiance and coherence statistics. Plane wave propagation differs from 
spherical wave propagation, as seen in the Rytov-Tatarski theory (Tatarski, 
1961), while beam wave propagation (Gaussian intensity profile) exhibits an 
intermediate behavior (ishimaru, 1978). The source must also be chosen to 
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keep the lateral extent of the source propagation within tfie spatial and spatial 
frequency limits of the computer simulation. 

Coherent plane, beam, and spherical waves differ In the shape of the 
isophase surfaces of the E-field. Beam wave sources possess a quadratic 
phase (tshimaru, 1978), 

Mp) ■ p*. <»•) 

where Rg represents the radius of curvature of the wavefront. Plane wave 
sources have an infinite radius of ojrvature so that 

♦(p) » oonstant. C**®®) 

Spherical wave sources have a radius of curvature equal to the propagation 
distance from the origin, or focus. S 

♦(p) ~ ^ (1®^) 

Turbulence introduces random phase shifts across the wavefront while 
diffraction and interference further distort the wavefront during propagation 
causing the light to become partially coherent. The resulting partially coherent 
E-field no longer possesses a simple isophase surface and the coherent 
wavefront characterizations no longer apply. The differences between plane, 
beam, and spherical propagation appear in differences in statistical properties. 
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such as normalized irradiance variance (Tatarski, 1961) 


^2 ^J23C^ pkne 

\).487 Cj spherical 


( 102 ) 


(see Ishimaru, 1978 for the more complicated beam wave expression), and the 
coherence length (Fried, 1966, p. 1380-1384) 

3.02 (Cj)"^ k-^ spherical 

ro = { (^03) 

1.68 k^ plane. 

When the E-fieid has propagated through turbulence to the far field of 
the source, diffraction has caused the E-field to expand laterally so that its 
statistical properties approach those of a spherical wave. The far field begins 
at a distance (Saleh and Teich, 1991) 

^ -10 ( 10 «) 


where p is the radius of the source. Spherical propagation properties result 
when the maximum radius of the source approximately satisfies 

For a stratospheric propagation scenario with L~200 km and \ = 500 nm, the 
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source radius must be less than 10 cm to obtain enough divergence to prochjce 
spherical propagation statistics at z - L. 

Plane wave propagation statistical properties result when the E-field is 
evaluated in the near field of the source because the E-field does not have the 
opportunity to diffract or expand significantly. Starting with Eq. (105) for the far 
field point, near field propagation satisfies 

z.±sL« 10^, (10S) 

10 X X 


or. 

> ^id X z. (“*07) 

For the stratospheric propagation scenario, source radius > 1 m for near 
field propagation. However, validity of the Fresnel approximation to the 
Huygens-Fresnel equation places an upper bound on (Saleh and Teich, 
1991) 

(p 5«)“ « 4 X, (lO®) 

(derived from considering the Taylor series expansion of the Cartesian 
expression for r). For the L = 200 km, X = 500 nm stratospheric propagation 
scenario, the upper bound on source lateral extent becomes < 350 m. A 
1024x1024 grid with the grid element size Ax = 0.99 cm from Eq. (92) gives a 
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maximum grid radius of only 5 m, so tiiat the Fresnel approximation is satisfied 
for any source represented on this grid. 

Ishimaru (1978) evaluated the irradiance statistics for the intermediate 
Gaussian beam wave case. The normalized irradiance variance along the 
beam axis depends upon the beam waist size, Wg, and the radius of curvature, 
Rq. Numerical integration of the log amplitude variance formula (Walters, 1994) 
provides a smooth transition from spherical wave variance statistics to plane 
wave statistics, with a dip in between, as shown in Fig. (10). Numerical 
simulation results with Gaussian and Airy-type sources of varying widths have a 
corresponding behavior, as shown in Fig. (11). 

The finite grid in a computer simulation places limitations on the source 
E-field. The physical grid element size sets a lower bound on the width of a 
narrow source approximating a point source. A source of width dose to a 
single grid element may still be undersampled, causing the resulting propagated 
field to exhibit sidelobes from absence of the proper high spatial frequences in 
the representation. Correspondingly, a spherical wave from a very narrow 
source will propagate with large angular divergence that may exceed the 
physical dimensions of the grid and cause energy to leak off the grid and be 
aliased back into the field. Thus, the source must be wide enough to constrain 
the propagated field so that most (>90%) of the energy remains within the grid. 
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Normalizad irrodionc* vorionee / Bcto squorad (tpharieai) 



Q 


Figura 10 Normalized irradiance variance versus beam waist size Wg, wl 
Q^icWg’/XL, from numerical integration of Ishimaru's predictions for beam 
waves. 
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Normalized irrodionce vorionce / Beta squared (spherical) 



Figura 11 Normalized irradiance variance versus beam waist size Wq, where 
rbffWo^/XL. Solid line = Airy-type source; dashed line = Gaussian source. 
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— >- .«UM 11. 'JUKI, A.^^»|*«JU|.NiJ^.|||UWslW!H*WJ 

Similwty. a plane wave source cwmot extend too near the grid edge because 
turbulence can scatter energy off the grid only to be aliased back in. 

Statistical calculations on the final beam pattern require a reasonably 
uniform central patch to perform computations. If the region over which 
statistical calculations are made has variations, such as sidelobes arising simply 
from vacuum propagation, these variations become a part of the statistics, such 
as the normalized irradianoe variance. Often, Gaussian type profiles present a 
minimum of such variations because they do not have as much energy in high 
spatial frequencies. However, because they are rK>t fiat over the calculation 
region, different regions of the beam must often be weighted with respect to 
their mean intensity in doing statisticai calculations such as the normalized 
irradiance variance. 

To meet these constraints, Martin and Flatt4 (1990) applied a quadratic 
phase curvature to a Gaussian source, thereby increasing the divergence and 
flattening the final irradiance pattern. Spherical wave sources were simply 
chosen narrower than the plane wave sources to make them diverge more. 
Flatt6, Wang, and Martin (1993) also used a super-Gaussian to model an 
extended beam source. 

The simulations presented here used an alternate method suggested by 
Ellerbroek (1993). The final field at z - L was specified as an aperture of 
radius equal to one half the grid radius and given a quadratic phase 
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corresponding to a spherical wave originating at the origin z = 0 (Fig. (12)). 

This aperture field was then backward propagated without turbulence from z = L 
to the focus z > 0, effectively Fourier transforming the field and yielding Vne 
numerical equivalent of the Airy pattern (Fig. (13)). This source was then used 
as the approximation to a point source for all spherical wave propagations. The 
advantages of this method were that the source was quite narrow, having 
appreciable amplitude over only a few grid elements, and that the zero 
turbulence propagated field at z - L was necessarily constrained within the grid 
and possessed a central region with uniform illumination. 

The significant energy at high spatial frequencies required to represent 
the sharp edges of the initial aperture at z - L created one drawback since it 
caused strong Fresnel fluctuations at intermediate distances between 0 and L. 
To investigate the significance of these Fresnel fluctuations, the energy at these 
high spatial frequencies were reduced by windowing the aperture before the 
backward propagation by applying a Gaussian rolloff to the cylinder edges. 


1 , r < 0.7 r' 

I . ri0.7r', n=2, 
0.3 r J ' 


(109) 


where r* represented the corresponding aperture radius. Because this rolloff 
reduced the energy in high spatial frequencies for the source representation, 
this slightly increased the size of the final central region of the E-fietd where 
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Figure 12 Cross section of intensity profile for final propagated beam showing 
confinement of energy and uniform central illumination (beam radius = 70 grid 
elements here). 



Figure 13 Cross section of intensity profile of numerical equivalent to Airy 
source. 
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statistical calculations could be performed. Other values for the power n on the 
exponent were tried (n s 1, n s 5 , n = 10) but the simple Gaussian (n = 2) 
worked best. While this rolloff appeared useful, identical runs, one with the 
sharply defined aperture edge and one with a rolled-off Gaussian edge, showed 
less than 0.5% difference in the normalized irradiance variance values. This 
difference was negligible compared to statistical fluctuations in normalized 
irradiance variance of up to ^10% between runs ^th different random phase 
screens. Empirical observation of simulations revealed that once even modest 
amounts of turbulence existed (typically ^ within the Rytov 

regime), the high spatial frequency energy introduced by the turbulence 
dominated the details of the source. 

Another deficiency was the representation of the initial cylindrical 
aperture on a rectangular grid. Due to the Cartesian nature of the grid, the 
curved beam edge was actually jagged instead of circular. However, for large 
enough grid (e.g. 1024x1024) this departure from a cylinder introduced 
negligible effect. The only case where a difference was noted was with a 
256x256 grid, and increasing the beam radius a small fraction of a grid element 
size eliminated that difference. 

Another limitation of this method was that the final E-field could not 
approach the grid radius. If this happened, then high turbulence strengths 
would scatter energy off the grid, producing aliasing. Due to the periodic 
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Fourier transform method of propagation, this energy would not be lost, but 
aliased back in at lower spatial frequencies. Martin and Fiattd (1990) 
implemented an attenuating region just inside the grid radius to absorb this 
energy and to prevent its being aliased back in. This obviously introduced a 
type of windowing function in the spatial domain and did not conserve energy in 
the field, possibly complicating final field comparisons. Rather than include 
those effects, these simulations simply chose the aperture radius suffidentiy 
small (at one half the grid radius) to prevent significant energy scattering off the 
grid in the spatial domain while still providing a relatively uniformly illuminated 
central region for statistical calculations. 

Obviously, other choices for a source existed. For example, the initial 
field could have been expressed according to an analytic expression for the 
Airy pattern (Fig. (14)). This source gave a vacuum propagated final field that 
was very dose to a broad cylinder, but which exhibited noticeable sidelobes at 
the edges (Fig. (15)). These sidelobes came from spatial truncation of the Airy 
pattern. The back propagated numerical Airy pattern of Fig. (13) differs from an 
analytic Airy pattern in a way that eliminates the sidelobes in the final irradiance 
field. 

A point source could be modeled by a single, nonzero pixel at the center 
of the grid being nonzero, as shown in Fig. (16). Figure (17) shows that the 
vacuum propagated field exhibited noticeable sidelobes at the edges and 
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covered the entire aperture, meaning that turbulence would immediately scatter 
energy out of the grid only to be aliased back in and distort the field. A similar 
source was a uniformly illuminated grid at z > L with quadratic phase, backward 
propagated to z = 0. While this eliminated the vacuum propagation ringing 
problem by definition, it immediately suffered from aliasing when turbulence was 
nonzero and scattered energy off the grid. 

E. MAXIMUM TURBULENCE STRENGTH 

One of the more difficult choices in a computer simulation of propagation 
is determining the range of turbulent strengths over which a simulation is 
valid, in general, the processes of optical propagation through a turbulent 
medium are not band limited in spatial frequency. But when a computer 
simulation uses a finite grid to implement a source function, to propagate via 
the Huygens-Fresnel prindpie, and to model the turbulence by phase screens, 
aliasing inevitably occurs due to the finite sampling interval. This problem only 
becomes more severe as turbulence strength increases. The coherence length 
measures the physical distance over which the mutual coherence of the E-field, 

(E* E ), dedines to e ’ of its peak value. As turbulence increases, the E-field 
fluctuates significantly over smaller and smaller distances and the coherence 
length decreases. An NxN discrete grid representation of the E-field samples at 
a specific minimum distance, and hence at a maximum spatial frequency. If the 
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E-field ^uctuates in less than this minimum distance, the grid samples will not 
represent the E-field accurately, causing aliasing. The question is not whether 
aliasing occurs, but how much aliasing occurs for a given set of simulation 
parameters, and at what point does aliasing invalidate the results of the 
simulation. 

These investigations identified three telltale signs of aliasing. First, as 
aliasing became significant, irradiance plots of the turbulent E-field lost structure 
and exhibited an isotropic fine-grained appearance, as Fig. (18) shows. 

Second, the irradiance that should have been roughly radially symmetric 


Figure 18 Intensity plot showing fine-grained pattern due to significant aliasing. 
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became bounded within a fuzzy rectangular region, as Fig. (19) illustrates. The 
irradiance at the center of the image increased, or started peaking, while the 
irradiance at the edges decreased, as shown in Fig. (20). 

Figure (21) also illustrates this peaking behavior by plotting average 
irradiance versus radius, in units of grid element size. The E-field propagated 
with zero turbulence had the radially symmetric aperture profile in this 
simulation. As the turbulence became stronger, the sharp edges of the initial 
cylindrical irradiance pattern became rounded and energy spread outward. 
However, when significant aliasing started occurring around ~ 5 for this 
64x64 ghd representation, the energy began creeping inward and the center of 
the field started increasing in irradiance. 

The irradiance at the edge of the grid was -1/40 that of the center when 
significant aliasing began. This indicates that very little energy leaked off the 
spatial grid due to beam divergence. Actually, the use of Fourier transforms in 
the propagation preserved the energy so that any energy that leaked off the 
grid was aliased back into the field on the grid. Energy that leaks off in the 
spatial domain results from undersampiing in the spatial frequency domain, and 
vice versa. 

Aliasing most often occurs from undersampiing in the spatial domain, 
which corresponds to energy leaking off the grid in the spatial frequency 
domain. Prudent choice of source and the final irradiance patterns across the 
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grid can minimize the amount of energy tfiat leaks off the grid to be aliased 
back in for the spatial domain. However, as turbulence strength increases, the 
coherence length decreases and more energy occurs at higher and higher 
spatial frequencies. Eventually spatial undersampling occurs and this high 
spatial frequency energy leaks off the spatial frequency grid and is aliased back 
in. Figure (22) gives the radial power spectral density corresponding to Fig. 

(21) and shows the spread of energy to higher spatial wavenumbers as 
turbulence strength increases. The low turbulence spectra possess a strong 
central peak that becomes more rounded and flatter as turbulence strength 
increases, causing more energy to leak off the grid. Even* jally, enough energy 
has leaked off the grid and aliased bade in to make the spectrum approximately 
uniform at all spatial frequencies. 

Figure (23) actually shows this energy being reflected back in after it 
spills off. To reveal this phenomenon, the radial power spectrum for a 512x512 
grid has been divided by the corresponding part of the power spectrum on a 
1024x1024 grid. Since the larger grid has been chosen with a finer mesh than 
the 512x512 grid, the larger grid will not experience significant aliasing as soon. 
As turbulence strength increases, the ratio shows the enhancement of energy 
at the edge of the 512x512 grid (bins 200-256) as the energy leaks off / reflects 
back in. At very high levels of turbulence, the aliased energy has spread 
inward over all spatial frequencies. 
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Power 



FIgm 22 Power spectrum radial profile as turbulence strength increases: p, 
- 0 (solid line), 0.5 (pluses), 1.5 (asterisks), 5 (dots), 15 (diamonds), 50 
(triangles), 150 (squares), 500 (X's). 






0 so 100 ^50 200 250 

Radiol wavenumber (bin |) 


Figure 23 Power spectrum radial profile for 512x512 grid, divided by 
1024x1024 profile, showing energy aliasing back in. 0.15(plus), 
0.5(asterisk), 1.5(diamond),5(triangle). 15(square),50(X),150(dot). 
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Aliasing sppsars to causs tha intansity pa^ng bahavior of tha cantrai 
fiaid. Thase invastigations paramatarized this onset of peaking with five 
methods, which provided guidelines for maximum turbulence strength 
valid for a given grid size. An irradiance ratio and a Strehl ratio were first used 
to identify the onset of the peaking. A Gaussian source was propagated at 
turbulence strengths ptt* = [5x10^. 5x10-®. 5x10-*, 0.15, 0.5, 1.5, 5, 15, 50, 
150, 500] using 64x64, 128x128, 256x256, 512x512, and 1024x1024 grids. 

Ten runs for each case were used to calculate the average irradiance as a 
function of radius. These irradiance profiles were used to calculate an 
irradiance ratio, defined as the irradiance at the grid center divided by the 
irradiance at the maximum grid radius, and also the Strehl ratio, defined as the 
irradiance at the grid center for the E-field propagated through turbulence 
divided by the center irradiance for an E-fieid propagated through zero 
turbulence. Figure (24) plots the irradiance ratio versus turbulence strength for 
different grid sizes, and Fig. (25) shows the Strehl ratio versus turbulence 
strength for different grid sizes. 

As turbulence strength increased, diffraction and scattering spread the 
energy outward and caused the irradiance and Strehl ratios to decrease. 
Eventually, both ratios reached a minimum and then started increasing because 
of the alias-induced intensity peaking. This minimum occurred at higher 
turbulence strengths for larger grid size because the larger grid sizes sampled 
with a finer mesh and did not significantly alias as soon. Using estimates of 
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Figura 24 Irradiance ratio versus turbulence strength ( = [5x1 5x10^ 

for grid sizes; 64x64 (asterisks), 128x128 (diamonds), 256x256 (triangles), 
512x512 (squares), 1024x1024 (X's). 
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these minima as a guide to the onset of significant aliasing, Fig. (26) plots 
these Po’ values versus grid size from the irradiance and Strehl ratios, along 
with least squares fits and extrapoiations to larger grid sizes. 

Martin and Flatty (1988) and others have validated computer simulations 
for predicting statistical properties such as normalized irradiance variance from 
the E-fieid propagated through turbulence. The departure of the computer 
simulation E-fields and intensities from their known or expected smooth 
behavior as turbulence increases represents another telltale sign of significant 
aliasing. Figure (27) shows the normalized irradiance variance calculated from 
computer simulated E-fields versus using grids of sizes 64x64, 128x128, 
256x256, 512x512, and 1024x1024. The irradiance value used for 
normalization was taken as the average irradiance over the central calculation 
region. Though this single value normalization is not optimal (see the following 
section, F. AddWonal Shnulaflon Parameters), it does make the normalized 
irradiance variance calculation sensitive to the peaking behavior because 
peaking adds large amounts of unphysical variance. In the Rytov regime 
(3o^ ^ 1.0), all grids successfully simulated the E-field as indicated by their 
normalized irradiance variances agreeing with the Rytov-Tatarski theory. In the 
saturation regime (Pg^ ^ 1.0) the normalized irradiance variance saturates and 
turns downward, eventually approaching unity according to asymptotic theory. 
Assuming that the 1024x1024 grid provides the most accurate propagation 
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Figure 27 Normalized inradiance variance versus turbulence strength from 
simulations with different grid sizes. 
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simulation, the 64x64 grid normalized irradiance variance departs from the 
1024x1024 values before the peak was reached, and 128x128 grid normalized 
irradiance variance departed just beyond the peak. Their meshes were not 
small enough to adequately sample the E-field, and the resulting aliasing 
caused peaking that drove up the normalized irradiance variance calculation. 
The 256x256 grid produced the saturation peak, but soon suffered from 
significant aliasing. The 512x512 calculations successfully produced the peak 
but departed from the 1024x1024 predictions before Po* ” 50 . Judging from 
these behaviors, the 1024x1024 grid probably remains valid through 3,,^ > 50 . 
All grids eventually showed the normalized irradiance variance anomalously 
rising for high enough turbulence strength because aliasing produced peaking. 
Estimates of the maximum values at which the five grid sizes remained valid 
are plotted, along with a least squares fit extrapolated to larger grid sizes, in 
Fig. (28). 

In a similar manner, the coherence length and the half width at half 
maximum (HWHM) of the atmospheric MTF predicted using simulations with 
different grid sizes also show unphysical behavior as aliasing became 
significant. Figure (29) plots the simulation derived coherence length, 
normalized by the Kolmogorov turbulence theoretical values, versus turbulence 
strength measured by p^^. Figure (30) shows the corresponding plot for 
HWHM. Both show the eventual strong rise in coherence length and HWHM 


76 


Maximum Beto squared 







0.001 0.010 0.100 1.000 10.000 100.000 1000.000 
Beto squored 


Figure 29 Spherical wave cohefence length from simulations versus turbulence 
strength for grid sizes 128x128 (diamonds), 256x256 (triangles), 512x512 
(squares), and 1024x1024 (X's). 
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Figure 30 Spherical wave HWHM from simulations versus turbulence strength 
for grid sizes 128x128 (diamonds). 256x256 (triangles). 512x512 (squares), and 
1024x1024 (X’s). 
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compared to theory for strong turbulence when aliasing produces peaking. 
Estimates of the minima of the plots can be used as another indicator of the 
turbulence strength at which significant aliasing occurs. Figure (31) plots these 
points versus grid size, along with a least squares estimate extended to larger 
grid sizes. 

The onset of peaking and departure from expected physical behavior 
provide symptoms of aliasing but do not indicate the amount of aliasing 
required to cause them. Determination of the fraction of total field energy that 
is aliased serves as one parameterization of the amount of aliasing and 
accomplishes two goals: (1) relate the amount of aliasing occurring to an easiiy 
measurable characteristic of the simulation , and (2) determine how much 
aliasing must occur to invalidate ttie computer simulation. 

First, to parameterize the amount of energy aliased, the same size 
Gaussian source was applied to 64;^, 128x128, 256x256, and 1024x1024 
grids and then propagated through turbulence. The width of the source was 
chosen so that the Fourier transform had approximately the same width on the 
final grid and the Gaussian sources on the different grids were normalized to 
the same energies. Since identical sources were applied in the spatial domain 
according to the analytical Gaussian formula, the spatial frequency 
representation of these sources improved as the grid size increased and the 
mesh became finer. The 1024x1024 grid served as an approximation to infinite 


80 





Moximum Beta squared 



10 100 1000 10000 

Grid Size 


Figure 31 Estimates and least squares fits from spherical wave coherence 
length (triangles / asterisks) and HWHM (diamonds / pluses) of the maximum 
turbulence strength versus grid size for valid simulations. 





grid size. Due to its larger extent and finer mesh, the 1024x1024 spatial 
frequency representation of the Gaussians contained spectral energy in the 
region outside the spectral footprints of the smaller grids, as Fig. (32) illustrates. 
Since each grid started with the same total energy, the fraction of spectral 
energy in the 1024x1024 grid lying outside the spectra) footprints of the smaller 
grids approximated the amount of energy aliased in the smaller grids. This 
process was then applied to propagations of the Gaussian beams through 
increasing strengths of turbulence, well into the significant aliasing regions for 



Figtae 32 Illustration of the spectral representations of the Gaussian beam 
with different grid sizes. 
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the smaller grids. This method provided estimates of the fraction of energy 
aliased versus turbulence strength for the three smaller grids. 

To achieve the first goal of relating the amount of aliasing to a 
measurable simulation parameter, the average radial power spectral density of 
the propagated E-field was calculated for each turbulence strength with the 
smaller grids. Using these profiles, a spectral ratio was calculated, defined as 
the ratio of the power spectral density at the grid center to the power spectral 
density at the maximum grid radius. Figure (33) plots these results, and also 
includes spectral ratios for 512x512 and 1024x1024 grids. Using the 
information on fraction of energy aliased versus and spectral ratio versus 
Po^. the turbulence strength p,,^ served as the connection between fraction of 
energy aliased and the measurable simulation parameter of spectral ratio. 

Figure (34) shows plots of these spectral ratios versus corresponding fraction of 
energy aliased. The log-log plot behavior proved approximately linear over the 
range 0.001 to 0.1 , which appears very useful because 0.1% of energy aliased 
probably does not affect the simulation results significantly while more than 
10% of energy aliased probably does. 

Using the data points of Fig. (34) and performing a linear least squares 
fit to the logarithms 

loatoCfl)' eiog«l«5) * ioBioA (ii«) 
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or equivalently, 


/? - 10^ {FRi^, O'***) 

where R refers to spectral ratio and FR refers to fraction of energy aliased, 
gives A- -1.1 ± 0.1, B- -2.1 ± 0.1 . Substituting the A and B values into Eq. 
(Ill) and rearranging. 




1 

yiwnii 


( 112 ) 


This equation achieves the first goal of relating the fraction of energy aliased to 
a measurable quantity from the simulation. 

However, Eq. (112) should be used with caution because it appears to 
be specific to the narrow Gaussian source used to derive it. Figure (35) shows 
the same calculations carried out for a source that was the Fourier transform of 
an aperture of radius equal to 3/4 of the 64x64 grid radius (i.e. Airy-type 
source). The linear region did not appear, though the plot matched the 
Gaussian source plot when the fraction of energy aliased was greater than 
0.1 . The differences arise from the different spectral representations of the 
propagated E-fields. The abrupt, cylinder-shaped irradiance patterns start with 
more energy at higher spatial frequencies in the spectral domain than the 
Gaussian patterns and thus never get below 0.001 fraction of energy aliased 
even at very low turbulence strengths and high power spectral ratios. 
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Flgim 35 Power spectral density ratio versus fraction of energy aliased for 
Airy-type source; grid sizes; 64x64 (dotted line), 128x128 (dashed line), 
256x256 (solid line). 
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The sea>nd goal of determining the amount of aliasing that invalidates 
the simulation used the fraction of spectral energy aliased to generate another 
guideline for maximum strength of turbulence for a given grid size where the 
maximum fraction of energy allowed to be aliased in a simulation was 10%. 
Equation (112) predicts a spectral ratio of approximately 10 for 10% of energy 
aliased. Linearly interpolation between the data points of Fig. (33) provided 
estimates of the turbulence strengfris C„^ (thus 3^^) that gave a spectral ratio 
= 10, for 10% energy aliased, for the five grid sizes. Figure (36) plots these 
values and a least squares fit extended to larger grid sizes. The least squares 
fit as a function of grid size N was 

•ogioPo = 0.9 log,o( V) - 0.9. ('•‘*3) 

or, 

p5 = 0.1 

The fraction of energy aliased also provided estimates of maximum 
for the Airy-type source. Figure (37) plots the fraction of energy aliased versus 
turbulence strength for the 64x64, 128x128, and 256x256 grids. Again 
choosing 10% energy aliased and using linear interpolation yields data points 
for the solid line in Fig. (36) that shows the maximum turbulence strengths 
obtained from Fig. (37). For 10% energy aliased, the predictions for maximum 
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Figure 36 Least squares fit (pluses, clotted line) of maximum beta squared 
versus grid size for a Gaussian source that had 10% energy aliased and for an 
Airy-type source (asterisks, solid line). 
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for the three smeller grids using the Airy-type source agree within ~'20% with 
the predictions fror«i the Gaussian source. Lower values of frac^onal energy 
aliased do not provide such agreement. 

Another measure of maximum turtHjIence strength involved the 
coherence length and grid size. Aliasing occurs because the E-fieid fluctuates 
significantly over scale sizes smaller than the grid element size. Intuitively, 
some conriection should exist between coherence ler^ith of the E-field and foe 
onset of significant aliasfog. Fried's coherence length r, represents the distance 
over which the atmospheric MTF falls to exp(-3.44) » 0.032 (Fried, 1966, p. 
1360-1383). Equation (38), which relates to C.^ and Eq. (66), which 
relates r, to C.^ arvi Eq (92), which relates the element size Ax to N, 
when combined, /ield foe turbulence correspondir>g to a given grid size 

for a specific number (y) of Ax's per r, 

.O.m ,9ptm1ot/m¥9 

pj - I 

' 0429 .pivte m¥9 

Figure (38) plote these p,’ for y « 2. 2.S. and 3 for foe spherical wave case 
The 10% allaaed energy line, indicated by plus symbols, lies newest the 2.5 Ax 
per r, Mne The E-fieW must be spatteiy s«T^)led with 2.5 Ax per r, to limit foe 
firacbon of energy aNased to 10%, analogous to foe Nyquist criterion two 
samplM per cycle 
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(39) provides a combined plot of tt>e maximum turbulence 
strength that produced valid E-fields for a given grid size from all of the 
fore^ing measures isity ratio. Strehl ratio, normalized irradiance variance, 
coherence length, HWriM, fraction of energy aliased, and coherence length 
(y s 2.5). Most estimates tie witNn a factor of 2 of each other, and the 10% 
energy aliased line (dotted/pluses) given by Eq. (114) represents an 
approximate lower bound to those estimates based upon onset of significant 
aliasing. 

Figure (39) provided three s«,|r^. .4 corKlusions: (1) the computer 
simulations remained v^id until approximatei} 10% of the energy became 
aliased (achieving the ffrst goal), (2) using the 10% energy attased Kne, the 
maximum turbuierKe strength p,’ for valid E-flelds for s grid size M was 

■ 0 -* 

and (3) maxtonum valid turbuierKe sfrength corresporxJed to approximately 
2.5 grid elements per r, (based on prox^ity of the y - 2.5 line to the 10% 
energy aliased Hne). CocKluslon (2) implies that doubto>g the grid size N (with 
grid element size Ax given by Eq. (92)) sUghtty less than doubles the maximum 
Pq’ that the igrki can simirfate, and that these investigations, which use a 
1024x1024 grid, shO(4d be vaM up to p,’^ 50. 
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F. ADDITIONAL SIMULATION PARAMETERS 

Additional aspects of a computer simulation were addressed to ensure 
validity, including the number of phase screens, the method of normalizing the 
irrat^ance variance, the number of realizatkms (propagations through 
turbulence) required for representative statistics, arKi the width of the final 
irra<fiance field on the grid. 

The number of phase screens must be large enough to represent the 
turbulence accurately along the path and produce proper irradiance and 
coherence statistics. Martin and FlatM (1988) determ^ied the number of phase 
screens by requiririg that the variance due to propagation over the distance Az 
between phase screens be less than 1/10 the vwianoe from the total 
propagation over the distance L 

ofCAz) < ai oJ(L), (1^7) 

and additionally that the vMue of the variance from one step be less than 0.1 

o^{Ai) < 0.1 . 

WHh these consideratlorui. they generally used 20 phase screeiis for the^ 
simulations 

These investigations examir>ed the irradianoe and coherence statistics 
(firecfiy arKi concluded that approximately 30 phase screens were recMr*0 to 
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ensure simulation validity. Figure (40) shows a set of spherical wave 
normalized irradiance variance simulations for a 512x512 grid, ~ 1.5, and 
the number of phase screens varying between 2, 4, 8, 16, 32, and 64. 

Assuming the 64 phase screen case most doseiy approximated physical reality, 
as UfH as 8 phase screens gave a normalized irradiance variance within ~10% of 
this value. By 32 phase screens, the ncxmalized irradiance variance had 
stabilized to within 2% of the 64 phase screen value. For this level of 
titfbulence that lies at the beginnirtg of the saturation regime, the addition of 
phase screens beyond 32 did not affect the normalized irradiance variance 
significantly. 

Figure (41) shows the normalized ^radiance variance from propagations 
with 2, 4, 8, 16, 32, and 64 phase screens at 50 with a 1024x1024 grid 
and 30 realizations for each mxnber of phase screens. This strength of 
turbulence represents the limit of simulation validity with the 1024x1024 grid 
and lies in the strong saturation regime beyond the normalized irradiance 
variance peak. In this case, 32 phase screens provided a normalized 
irradiance variance within 3% of the 64 phase sc r ee n case, arKl the 16 phase 
screen value was within approximately 10%. Because 32 phase screens 
appears to increase the accuracy of the simutation by ~5% compared to 16 
phase screens, me larger number of phase screens was chosen for these 
investigabons 
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W t lji|. II 



Using th« sssne 1024x1024 resUzations, the coherence length and 
HWHM of the atmospheric MTF were calculated for the 2, 4, 8, 16, 32, and 64 
phase screen cases. Figure (42) plots the results. For coherence length, the 
32 phase screen value was within approximately 1% of the 64 phase screen 
value, and the 16 phase screen value was within approximately 3%. The 
HWHM plot indicates 2% ar>d 3% agreements for 32 and 16 phase screens, 
respectively. The 32 phase screens used for these investigations thus proved 
sufncient for cohererKe length and HWHM calculations. 

The method of normaltzing the irra<tonce variance and the number of 
realizations to use for statistical accuracy proved to be interrelated 
considerations Turbulence diffrects and scatters energy outward from an 
initialty well-defined beam For the A^-type source used in these 
investigations, the average irradiance over the central portion of the final 
propagated field was uniform for zero turbulence but became a function of 
radial cNstance from the propagation axis when turbi^erKe was present. This 
was an artifact of foe computer simulation that had to use a beam spatially 
cortfined to the grid rather than a true spherical (or plane) wave Figure (43) 
shows irra(fian<» versus radial (fistanoe r for foe carttral 256x256 portion of 
1024x1(^4 (fod simulations, averaged over 30 reaUzations, and using a 
turbulence strength fi,’ « 50 The average radM irradlWKe vvies by 
iHPproximafoly 15% ovw’ the widfo of this calculation region This rarfial 


99 




of ph«M BCTMlIt 


W pnM« tcTMOf at turt>utooce strangth 0o * 50 




Averoge intensity 








variation of the average irradiance was removed from the rKKmaltzed irradiance 
variance calculations with a method smilar to that of Flattd, Wang, and Martin 
(1993). The calculation region was restricted to half the radius of the zero 
turbulence irradiance pattern, which was previously chosen as one half the grid 
radius. WHh the 1024x1024 grid, this calculation region equaled the largest 
radius circle inside the central 256x256 portion of the grid. This disk was 
further divided into concentric rings arKl the irradiance in each ring averaged 
over ail 30 realizations. These average ring irradiar)oes were used to normalize 
the irradiance variar>ce calculated from each field Smaller ring size (thus more 
rings) reduced the number of poirrts in the irradianoe average for each nng and 
required more realizations to ensure a sufficient nunber of points to yield a 
stable average irradiance 

The ring width had to be smaH enough to compensate for the radial 
variation in average irradianoe and the number of realtzations targe enough to 
provide enough points to yield representative average values To determine 
suitable ring width and number of realizations. 50 1024x1024 realizations with 
the Airy>type source at » 1 . 5 , 3 , 10 . arKt 20 were run and the normakzed 
irradiance variances caloMted usirtg ring widths of 128, 32. 8 , 4, 2, and 1 grid 
elements, dx, and with number of re^zations in the average equal to 1, 5, 10, 
15, . ,50 Figure (44) plots the rest^ for the 89 * s 10 set The pluses lir>e 
corresponds to 1 ring of width 128 dx (i e the whole calculation region) and a 
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Nofmoli<«d irrodionc« 







8 tabi« rvormaltzed irr»tianc8 value was achieved using about 20 realizations. 
However, a single hrtg does r)ot compensate for the radial variation in average 
irra<^ar)ce and can thus give erroneous normalized irradianoe values (for the 
same reason, it becomes sensitive to the pe^ng behavior from significant 
aliasing) The X line represents the division of the disk into 4 rings of width 32 
Ax and re€M>^ ‘‘25 realizations in the average to give a stable normalized 
irradiance vwiance. SmaHer ring widths (Knes with squares, diamonds, 
triangles, and asterisks lines) all showed similar behavior and required - 30 runs 
to reach a stable average normi^ed irradiance variance. The s 1.5, 3, 
arrd 20 runs aN provided similar results. A similar p,* s io series with a 
Gaussian source required a mir^mum of 8 rings (16 Ax each) and 30 
re^izations to active stable average normalized irradiance variarK^es. 
Conse<^ientty, these investigatiorrs used 30 reakzations as a guideline for valid 
simi4ation, and chose a rirrg width tow elements to irilow the greatest 
adjustment to real variations in average redial kiterrsity and yet remain 
computationally efficient. As Fig. (44) ^lows. insufficient averaj^ over enough 
realizations to yield truly representative average values can reduce the 
normalized irrackance varUmca by 15 • 30%. 

Sknkarty, stable coherence length values (discussed in section. H. 
Cobet e tice Lengffi) required averaging over multiple realizations. To determine 
the number of realizations requM fo the average, 50 realizations usir>g a 
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1024x1024 grid were g«n«r«t«d for p,’ s 3 and 20 and the central 256x256 
portion of the fields again used to calculate foe atmospheric MTFs. These 

individual MTFs were then averaged in groups of 1. 5, 10, 15.50 before 

calculating the coherence ler>gfos. Figure (45) shows foe resulting coherence 
lengths, r,. (These particular runs used L > 150 m and X = 420 nm, giving a 
FresTHil length ■ L / 1 % - 10 mm and resultir>g in millimeter-sized 
colwence lengths.) As few as 5 realizations in the average yielded coherence 
lengths within 4% of foe 50-faalization values. These investigations still used 
30 realizatiorM in the average because these 30 fields were already available 
from the normalized irradiance variance calcuiations. 

The irradiacKe ratio provided an indicator of the appropriate final beam 
radius to use to avoid excessive scatter of enwgy off the grid. The maximum 
Pq’ for a given grid varies with the final beam radius because larger radii scatter 
energy off the grid sooner and cause aliasing at lower To characterize this 
behavior, 64x64, 128x128, and 256x256 grid simulations were run at turbulence 
*<i’'4hgfos Pp' (5x10'*, 5x10^ for ffoal bewn radii of 4/8, 5/8, 6/8, and 7/8 grid 
radius, and the resulting irradiartce ratios (average irradiance at center divided 
by average irracttarKe at grid radius) were plotted versus ^ . Figure (46) 
shows the plot for the 256x256 grid. Again, the mir)ima correspond to foe onset 
of peaking and occur at lower Pp’ for larger final beam radius. Linear least 
squares fit of the final beam radii to these V values for the minimum 
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irradianoe ratios indicates that a final beam radius of approximately 0.7 grid 
radius corresponds with the 10% energy aliased cutoff for the grids. A final 
beam of 0.7 grid radius provided the largest illuminated central region while sfill 
meeting the 10% aliased energy criterion. To be somewhat conservative, these 
investigations used a final beam radius = 0.5 grid radius to reduce further the 
energy scattered off the grid. 

Q. PHASE SCREENS 

1. Pliaae Saaen OeneraSoii 

WHh the split-step method, the effects of turbulence along the 
optical path are introduced into the simulation by dividing the optical path into 
steps and applying a rarKfom phase to the complex E-fieid at each step. As 
long as the steps are smalt erKXjgh that geometrical optics approximately 
applies, the E-fleld only acquires a random phase change as it propagates 
across each step (Knepp, 1983). Diffraction as the field propagates across 
many steps then produces the ampfitude variations. The random phases are 
assumed to be Gaussian distributed about a zero mean with variance 
proportionai to the turbuienoe strength C,’ and possessing spatial structure 
function consistent with the assumed Kolmogorov turbulence and appropriate 
inner scale. Knepp (1963) and Martin and Flatti (1988) describe the process 
of generating the phase screen with these characteristics, as discussed below. 
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The phase screen generation begins in the spatial frequency 
domain by imposing the proper spatial structure function. An NxN grid of 


complex numbers 6 o(K„Ky) is formed whose real and imaginary parts are each 
Gaussian distributed random numbers with zero mean and unity standard 
deviation. This 6 o(K,.iCy) represents the Fourier transform of a grid of 
uncorrelated Gaussian distributed random numbers 6 o(x.y) representing phases. 
The proper spatial structure function corresponding to turbulence statistics is 
imposed upon the random phases 6 o(x.y) by applying a filter A(K,.iCy) to 

Taking the magnitude of both sides of Eq. (119), squaring, assuming that the 
filter function is real, and then takirtg expectation values gives 

< I e(V,Ky) I* > • ( I I* > • 

Where use was made of the fact that the Gaussian random numbers 6o(K„Ky) 

have a variance of 1. The two-dimensional power spectral density of the phase 

^s(k»S) rela^ to 0 (k„k,) by (Goodman, 1985) 

Where Ak represents the grid element size in the spatial frequertcy domain. 

The Hankei transform of the power spectra density F,(k) gives the phase 
structure function D,(p) that characterizes the spatial distribution of the phase 
fluctuations of the E-field (Tatarski, 1961) 
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0,M - /11 - 4(*p) 1 ''»(*.0) K d*. 

-«• _ 

where local iaotropy has been assumed. Tatarski derived the relation between 
the two-dimensional power spectral density of phase fluctuations F,(ic„Ky) wkI 
the tiwee-dimensional spectrum of index of refraction fluctuation-. <b,(K), 

- 2iilf* L ♦^(k ), (123) 

where x • ♦ kJ . This sequence (rf steps means that the phase screen 

can have the proper spatial statistics by starting with the proper spectrum of 
refractive index fluctuations (hence, the proper structure function). 

The spectrum of index of refraction fluctuations assuming 
Kolmogorov turbuienoe with irmer scale is 

- 0.033 C,*W «■”" P(ieW. (1W) 

Where F(k^ gives the inner scale dependence for kmer scale ^ (see Fig. 8). 
Sirt>stituting Eq. (124) into Eq. (123) gives the power spectral density of phase 
fluctuations 

F^{x^) - 2iclr* L (0.033) Ci{z) F(k^), (««) 

and using Eq. (121) spedfles the proper form for 6 (k,.k,), the corresponding 
Fourier transform of the random phases 

< I* > * (As)-* 2xk* L (0JK3) C!{J^ F(k^). (12«) 






Equation (120) than gives the corresponding filter function A(K.,Ky) to apply to 
the array of complex numbers 6 o(K„tCy) 

A(ic„Ky) » (Ax)-^ L (0.033) C!!{z) F(ic^). 

Since the Kolmof^rov spectrum (x has a singularity at k - 0, fiie 
point in the filter fur>ction is set to zero, removing overall piston (i.e. common 
phase offset over whole screen) from the phase (Cochrane, 1985) and keeping 
the spectral energy finite. Conversion to a discrete grid representation occurs 
by substituting x, = n. Ax, k, = n^ Ak. x* * (Ax)* (n,* ♦ n^*), arnl 
Ax 2ic/(N Ax), where Ax is the grid element size in the x-y domain and (n,.ny) 
are grid coordinates. The Fourier transform (FT) of the filtered array of random 
variables gives the phase screen in the spatial domain 

e(x,y) . 0.0864 Hr ^C^(z) L {N Ax)^ FT[ «,(/»„#»/ ]. 

Since this phase screen is actually complex-valued, both the real and imaginary 
parts represent vaNd random phase screens that were separately applied to the 
E-field (Cochrane. 1965). 

2. Low 8 pa6al Frequency Conecdon 

Due to the finite size of the the above phase screen will not 
have the proper structure function for separations of the order of the grid width. 
Low spatiai fiequency components, espedatiy tilt-type terms, are under- 
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represented (Cochrane, 1985). These computer simulations incorporate an 
algorithm formulated by Cochrane that employs an expansion of the phase 
screen in a Karhunen-Loeve basis set (whose components are mutually 
orthogonal) to correct some of the low spatial frequency terms. 

The general idea of the low frequency correction builds upon the 
above procedure to generate a phase screen. The low spatial frequency 
contribution to the phase screen is a superposition of orthogonal low spatial 
frequency terms, just like the superposition of complex exponential terms by 
Fourier transform in the phase screen generation process above. The strength 
of each low spatial frequer>cy term is a Gaussian random variable with an 
appropriate variance, just as the Gaussian raryJom numbers above were filtered 
to give the proper varUmce and thus determine the strength of the 
corresporKJing complex exponential in the spatial domain. Thus, the two 
objectives involve finding appropriate set of orthogor^al low spatial frequency 
functions, and determkiing the corresporKiing variances applinble to 
atmospheric turbulence. 

An arbitrary ftmction can be expanded in terms of Karhunen- 
Loeve hjnctions that are orthogonal by definition. To determkie a set of 
Karhunen-Loeve functions appropriate for Kolmogorov turbulence, Cochrane 
builds on the work of Non arKi (1975) considers expansion of an arbitrary 
function ^r,0) over a circular aperture of radius R in terms of Zemike 


112 






polynomialt (Bom and WoHd, 1970) 


♦(f.e) - Sa.z(p.e). 

whore p represents the normalized distance r/R, a, are the expansion 
ooefftcients, and Z, are the Zemike polynomials. With W(r/R) representing the 
aperture function, the coefficients we 

«/ - fd*' mOH) ♦(r,») Z,{rlRA). (1M) 

NoN assumed that these coefficients were (Gaussian random variables with zero 
mean arxl with a covariance 

(- /«!» fdf' tn())Wtp') ^(p.») < ♦(«?) ♦(wpO > 

Fourier transforming to the spatial frequency domain 

•//«<«' o;w •(ftiR.t'm («*) 

where Q^k) represents the Fourier transform of the |th Zemike polynomial, and 
^k/R,k'/R) represents the Kolmogorov spectrum of phase fluctuations. Non 
analyticaily performed the integrals to give a covariance matrix in which the 
terms represent the expected covariances due to Kolmogorov turtulerioe. 

Cochrane (1985) notes that the Zemike polynomials cannot be 
used to form an orthogonal e xpa ns io n of the turbulerKeHtetorted phase 
because the expansion coefficients are correlated, IrKficated by nonzero off- 
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diagonal elements In Noll's covariance matrix. However, the eigenvectors of 
the Zemike covariance matrix serve as a Karhunen-Loeve basis set. These 
eigenvectors can represent turbulence because they are not correlated, i.e. 
each eigenvector is formed by superposition of Zemike polynomials in such 
a way that the K, are orthogonal, satisfying the first objective. Additionally, the 
cofTesponding eigenvalue \ multiplied by (D/rj*^ (where D is the aperture 
diameter and r« is Fried's cohererKe lerHJth (Fried, 1965)) gives the appropriate 
variar>ce for ttuit K, spatial compor^ent corresponding with Kolmogorov 
turbulence, satisfying the second objective. 

Spedficalty, the low spatial frequerK:y contribution to the phase 
screen can be expanded in terms of these Karhuner>>Loeve components K, 
(Cochrane, 1985) 

Awb 

♦W • E T, 

P»1 

where represwits the rnimber of low spatial frec^jerKy KartHJnen^.oeve 
terms IrKluded and the coefticients are Gaussian random numbers with 
variarKa arn] s ca l e d by {Otrj^ to the specHtc strer^gth of turbulecKe used. 
The simuiatiorui use the first five terms (p^ « 5 ) as a compromise between 
completeness of tow spatial frequerKy corr e ction and oomputatkxMl efflderKry. 
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Figure (47) shows a surface plot of one realization of the first two terms of the 
correction, which are very close to x- and y-tilt (i.e. phase terms linear in x and 
y, respectively). Figure (48) shows a realization of the 3rd, 4th, and 5th 
correction terms, which resemble the wavefront aberrations associated with 
defocussing and astigmatism in a conventional imaging system. To implement 
the low spatial frequency phase correction, (1) the scalar product was formed 
between the initial Founer transform phase screen and each Karhunen-Loeve 
function K^, jiving the relative strength of that in the initial phase screen; (2) 
this amount of each spatial component was then subtracted from the phase 
screen; and (3) the component was then added back to the phase screen in 
the proper amount given by the product of a Gaussian random number with 
variance and the factor to scale to ttie particular strength of 
turbulence used. 

Cochrane's computer routines only calculate the Karhunen-Loeve 
correction terms over the largest ctrde that fits inside the calculation gnd, as 
shown in Figs. (47) and (48). Correction of the E-field over the entire computer 
simulation grid requires the Karhurien-Loeve terms be calculated over an area 
that is at least /? larger on each side than the field grid (Cochrane, 1985). If 
the E-field grid size N is chosen as a power of 2. then the Karhunen-Loeve grid 
must be 2Nx2N. These simulations however used only an NxN grid 
(1024x1024) because the Sun SparcStations could not handle the memory 
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Figm 47 Terms 1 and 2 of Karttunen-Loeve low spatial frequency correction 
to phase screen, represented in optical path length for wavelength » 500 nm. 
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requirements of the doubled grid (2048x2048) without a major revision of the 
propagation code. Since the source function was chosen to minimize energy 
scattered off the grid, very little energy fell in the uncorrected comers of the grid 
and the simulation remained valid. 

Noll's covariance matrix assumes Kolmogorov turbulence 
spectrum with zero inner scale making the Karhunen-Loeve correction terms 
strictly apply only to this spectrum. These simulations however incorporate 
nonzero inner scales into the spectrum of refractive index fluctuations. 

Because the first five Karhunen-Loeve terms cover scale sizes on the order of 
the grid widtft, which is much larger than the inner scale size the Karhunen- 
Loeve correction terms derived for zero inner scale remain valid for correcting 
noruero inner scale phase screens. 

Cochrane (1985) showed that such low spatial frequency 
correction greatly improved the phase structure function. Two terms corrected 
the structure function to within 10% of the theoretical Kolmogorov behavior, and 
five terms corrected to within 5%. compared with 30 -1000% discrepancies 
without any correction at low spatial frequencies. Figure (49) plots the 
normalized irradiance variances from computer simulations as a function of 
for zero Karhunen-Loeve low frequency correction, two-term correction, and 
five-term correction. Zero correction underestimated the irradiance variance in 
8)e Rytov regime by approximately 5%, and the two-term tilt-type correction also 
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Figura 49 Normalized irradiance variance, normalized by versus 
turbulence strength for zero (dotted), two^erm (dashed), and five-term (solid) 
Karhunen-Loeve low spatial frequency corrections to phase screens. 







underestimated by about 5%. The five-term correction with its focussing type 
terms raised the variance to within about 2 % of the theoretical variance. In the 
saturation regime, the two- and five-term corrections fit better. Figure (50) plots 
the coherence length rg, normalized by the theoretical coherence length, as a 
function of turbulence strength for zero, two-, and five-term Karhunen-Loeve low 
spatial frequency corrections. In the Rytov regime where simulation should 
dosely approximate theory, the zero correction overestimated the coherence 
length by approximately 35%. The tilt-type correction (two terms) estimated the 
coherence length within about 5%, and the five-term correction achieved 
agreement within about 2%. These behaviors formed the guideline that some 
type of low spatial frequency correction was required to achieve valid 
coherence lengths from computer simulation. 

This low spatial frequency correction method remained 
computationally feasible because the Noll covariance matrix only needed to be 
calculated once and because only a few terms of the Karhunen-Loeve 
expansion were used. However, the code can become memory intensive 
because it saves a full NxN grid for each Karhunen-Loeve function in addition 
to the NxN phase screen itself. Implementation with 2048x2048 or larger grids 
becomes problematical except on very large computers with about one gigabyte 
of RAM. Fried(1993) has proposed a simpler x- and y-tilt correction algorithm 
and indicates that this performs almost as well as including the higher order 
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dotted line » 0 inner scale, 0—term KL correction 
dashed line 0 inner seole. 2-term KL correction 
solid lino ■ 0 inner scole, 5-term KL correction 



Flgiim 60 Spherical wave coherence length, normalized by theoretical 
coherence length, versus turbulence strength for zero (dotted), two-term 
(dashed), and five-term (solid) Karhunen-Loeve tow spatial frequency 
corrections to phase screens. 
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Karhunen-Loeve terms, while being less memory intensive and faster than the 
Karhunen-Loeve correction method. Figs. (49) and (50) indicate that the tilt- 
only correction may systematically underestimate the normalized irradiance 
variance by about 5% and overestimate the coherence length of a spherically 
diverging wave by about 5% in the Rytov regime. Depending on the 
application, the Karhunen-Loeve correction to higher order terms may prove a 
useful refinement. 


H. COHERENCE LENGTH 

The coherence length of the E-field was calculated from the atmospheric 
MTF with the theory outlined in Chapter II. but the actual implementation of the 
MTF calculation and parameterization of the coherence length required 
careful consideration. The calculation methods that worked best for these 
investigations are discussed in this section. 

Equation (60) of Chapter II relates the atmospheric MTF and the 
spherical wave structure function to the coherence properties of the E-field 


MTF^ 


m) 

9 * 




( W{r') ) 


(134) 


Again, U(f) represents the E-field and W(/) represents the aperture function. 
The autocorrelations of U and W indicated in Eq. (134) could be done by 
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summing the complex products of the E-fields over multiple pairs of points or by 
implementing the autocorrelation via FFT techniques. Both versions were tried 
for these investigations and produced similar results, but the FFT version 
provided a much more thorough autocorrelation with greater computational 
efficiency. 

The FFT autocorrelation technique used the MCF introduced earlier. 


Equation (57) defined the mutual coherence function as 

MCF . ( 4/*(A.f,) >, (13S) 

which, for a single time = tj = t and assuming homogeneity, may be written 
as the spatial autocorrelation of the E-field 

MCFif") = f dt U*it) U(t + t»). ( 13 «) 

Substituting Eq. (136) into the Fourier transform identity, 

MCF{r') = / d^l / dt"MCF{t") (137) 

gives 

MCF{f') = f df I f <ff" I f dr (/•(/) U{M") I j (138) 

Rearranging ffie integrations, 

MCF{F) ^ f df if dr (/•(/) I / U{r*r") 11 (139) 

Changing variables to ^ = 7 ♦ 7", 7" » 8 - 7, dT" » 
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MCF(f') - f df if dt U*{t) I / a U{i) 1 (140) 

The inner integral is the Fourier transform of U(l) and is denoted by FTI U ]. 

The next innermost integral is the inverse Fourier transform of U*(/) and is 
denoted by IFT[ U* 1. Then, 

MCF{F) = f df FTI U] tFT\U*] (141) 

which is yet another inverse Fourier transform, symbolically written 

MCF{F) - IFT[ FT\U\ IFT\U*\ ). (**^2) 

Equation (142) expresses the autocorrelation (MCF) of the E-field in terms of 
Fourier transform techniques easily implemented with discrete Fourier 
transforms. 

This Fourier transform method of calculating the autocorrelation of a 
function is faster and more complete foan the more laborious technique of 
averaging the products of the E-field at point pairs. The product E*(/| )E(f 2 ) is 

complex, but because of the averaging that occurs in the autocorrelation, foe 
real part attains a stable value while foe imaginary part averages toward zero. 

For the FFT implementation on a 1024x1024 grid, the imaginary part is typically 
*10** while the real part ranges between 0 and numbers on the order of unity. 
The point pairs technique produces an equivalent real part but reduces foe 
imagiriary part down to only -'10'’. These investigations used foe FFT version for 
all autocorrelations. Additionally, the autocorrelations were normalized by their 
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zero lag values to ensure that the magnitude of the MTF,^ calculated by 
Eq. (134) lies between 0 and 1. 

Because the spherical wave strudure function refers to points across a 
spherical wavefront, U and W must similarly represent the E-field across a 
spherical surface. However, the computer simulation implements the E-field on 
a plane perpendicular to the axis of propagation and incorporates ttie spherical 
wave nature of the E-field by applying a quadratic phase curvature across the 
plane. To convert from this plane representation of the E-field in the simulation 
to a spherical surface representation appropriate for the autocorrelation 
calculations of Eq. (134), the quadratic phase factor across the plane must be 
removed from the E-field. This effectively assumes that the amplitude and the 
phase fluctuations of the E-field across the plane of the computational grid 
closely approximate the E-field across the spherical wave. This assumption is 
justified because the maximum physical separation of the true spherical 
reference surface from the plane surface of the grid is at most 1/50 the 
coherence length of the E-field. However, the removal of the quadratic phase 
curvature from the E-field of the grid pro>ms crucial to the autocorrelation of U 
because the quadratic phase ftictor undergoes ~130 multiples of 2n phase 
change between the center and the outer edge of the grid for these simulations 
and would otherwise completely obscure the actual E-field fluctuations. 
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Propagation with a plane or beam wave requires a different approach to 
remove the proper amount of phase curvature. For a plane wave, the 
wavefront coincides with the plane of Vhe grid so that no phase curvature needs 
to be removed. A pure beam wave (Gaussian profile) exhibits a spherical 
phase with a radius of curvature larger than the propagation distance L, and 
this phase could be calculated analytically and removed. However, the use of a 
finite E-field confined to the propagation grid introduces phase effects other 
than simple quadratic curvature. Fortunately, the E-field propagated through 
zero turbulence contains this phase curvature information (Walters, 1994). For 
spherical wave propagations, the removal of the phase curvature by the 
analytical calculation and by using the zero turbulence propagated field 
provided identical coherence lengths. These investigations used the latter 
method for the coherence length calculations for all beam-like and plane-wave- 
like propagations. 

Equation (134) requires that the autocorrelation of U must be averaged 
over multiple realizations to achieve the long exposure MTF. To implement this 
requirement, the autocorrelations from several realizations were calculated and 
averaged together, and a coherence length then determined via Eq. (134). This 
method produced coherence lengffis that agreed with theory within ->5% in the 
Rytov regime. Recalling Fig. (45) for ffie coherence lengff) versus number of 
realizations used in the average, the number of fields included in the MTF 
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average may be as few as 5, though 20 provided a more statistically 
reproducible coherence length. Again, these simulations used 20 realizations 
because these fields had already been generated for the normalized irradiance 
variance calculations. 

Figure (51) plots the coherence lengths calculated from the average MTF 
and also plots the average of coherence lengths calculated with single 
realization MTF's. The average of individual realization coherence lengths 
exceeded the averaged MTF coherence length by '>20% at low turbulence 
strengths but eventually agreed within 5% near the saturation regime. The 
single realization MTF's showed a larger coherence length at low turbulence 
strengths because the contributions from tow spatial frequency components had 
not been reduced through averaging. At these low turbulence strengths, 
multiple realizations were required to average these low spatial frequency 
contributions and to achieve the appropriate long exposure MTF. At higher 
turbulence strengths near saturation, the E-field had more energy at high spatial 
frequencies that dominated the MTF. The autocorrelation for a single 
realization now averaged over many coherence lengths and yielded an MTF 
dose to the average MTF for multiple realizations. 

Using a finite beam to approximate a spherically diverging wave 
produced a radial dependence of the average irradiance that affected the 
coherence length calculation, but only to a small degree. To compensate for 
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mis radial deparKlenca, each E-field was divided by the radial average E-field 
magnitude from the 30 realizations. Similar to the irradiance variance 
calculation, this average E-fieid was calculated by dividing the grid into rings 
one grid element wide, taking the square root of the average intensity for each 
ling over the 30 realizations, and then performing an area weighted, running 
mean across the rings to smooth the variations. While this radial compensation 
proved essential for the normalized irradiance variance calculation, it only 
changed the coherence iengm by ~1%, whicm was significantly less thun the ~5% 
discrepancy from the low spatial frequency correction. 

The above considerations allowed calculation of the right-hand side of 
Eq. (134) for the atmospheric MTF. Because of the statistical nature of the 
propagation through turbulence, no set of realizations yielded an atmospheric 
MTF that exactly followed the exponential rolloff with distance predicted by the 
left-hand side of Eq. (134). Methods had to be developed to parameterize 
these atmospheric MTF's from the simulations arid extract an appropriate 
coherence length oorrespondirig to the structure tijnction D(r). 

For the case of Kolmogorov turbulefice, Fried (1966, p. 1380-1383) 
derived the wave structure function arid expressed it in terms of the single 
cohererice length parameter r, 

Oir) > Ij-f, (143) 
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where n = 5/3. The atmospheric MTF then becomes 


iUTF, 



( 144 ) 


To extract the coherence length ro arxt the exponent n (allowing the possibility 
that n may vary), take the natural logarithm of both sides and rearrange 
(Walters. 1993) 

,14.) 

Taking the natural logarithm again. 

ln( - ^ kHUTF^ir ))) = n h(f) - n ln(r,). (14«» 

This has the linear form y s ax b, where y represents the left hand side, 
a = n, X > in(r), and b - -n In(ro). The atmospheric MTF can now be 
characterized with two parameters, r, and n, or just the single parameter rg 
assuming n ■ 5/3. 

To imfNement these parameterizations. the atmospheric MTF was first 
calculated from the E-field using the autocorrelation methods above and then 
radially averaged to yield MTF,g.^(r) for use in Eq. (146). A linear least 
squares fit calculation provided the slope a - n and the intercept n In(ro), giving 
rg. To obtain the single parameter characterization, n was set to 5/3 in 
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Eq. (146) and least squares techniques applied to obtain the intercept and thus 
the single parameter Tq. 

Figure (52) shows the coherence lengths calculated with these least 
squares methods, divided by the Rytov theory coherence length for Kolmogorov 
turbulence. The two parameter characterization with both exponent n and 
coherence length rg never provided a consistent, smoothly varying coherence 
length. At low turbulence strengths in the Rytov regime where the E-fieid 
coherence length is larger than the calculation aperture, the two parameter fit 
predicted coherence lengths up to 50% higher than the theoretical value and 
thus appears unreliable. Additionally, it showed an anomalous bump around 
- 0.5 . The single parameter least squares technique with n = 5/3 
accurately characterized the coherence length within 5% at low turbulerK» 
strength, but still showed the bump at ~ 0.5 . 

An alternate technique to obtain a single parameter characterization 
assumed n = 5/3 in Eq. (144) and used a binary-type search, or iterative fit, to 
find the coherence length rg that minimized the variance between the average 
MTF,g^(r) and the right-hand side of Eq. (144). Though not analytical, the 
resulting rg gave an MTF that often fit the actual MTF,g^(r) more closely by eye 
than the least squares methods, especially for low turbulence where coherence 
lengths were larger than the grid size. To implement the technique, an initial 
large range of rg (for example, 0 to 100 m) was divided in half, the midpoint rg's 
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Figura S2 Coherence lengths calculated by one-parameter (dashed line) and 
two-parameter (dotted line) least squares techniques, the iterative fit technique 
(solid line), and the HWHM (dash-dot line). 
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of each half weie substituted into Eq. (144). and the variances were calculated. 
Whichever rg provided the least variance became the middle of the next range, 
and the other rg became the new high or low boundary. The process was then 
repeated, so that each iteration reduced the range of possible rg by 1A3. This 
procedure was iterated 100 times or until the variance was less than 1x10*'. 

Figure (52) shows the coherence lengths predicted by this iterative fit 
method along with the coherence lengths from the one* and two-parameter 
least squares fits, all divided by the Rytov theory coherence length for 
Kolmogorov turbulence. The iterative fit coherence length was ^-5% too large at 
low turbulence but did not show the anomalous bump around Pg^ ~ 0.5 . 

In the saturation regime where turbulence is high and/or path lengths are 
long, multiple scattering becomes significant and saturates the irradiance 
variance (Martin and Flatty, 1988). Correspondingly, this physical phenomenon 
may also affect the coherence length and modify the value of n or the form of 
Eq. (144) for the saturation regime. The half width at half maximum (HWHM) 
provided a coarser, one parameter method of characterizing the atmospheric 
MTF that did not depend on any assumptions about the form of the structure 
function and could be calculated directly from the atmospheric MTF by linearly 
interpolating between the pair of points bounding MTF,^„„,(r) = 1/2. Figure (52) 
plots the HWHM, divided by the HWHM predicted by Rytov theory assuming 
Kolmogorov turbulence. The HWHM plot starts at Pg^ = 0.05 because lower 
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turbulence strengths gave a coherence targe enough that the MTF,^(r) had 
not reached half its maximum value within the calculation region. The HWHM 
lengths followed the theoretical values ^thin 5% in the Rytov regime and did 
not show the anomalous bump around ^ ~ 0.5 . Note that the iterative fit 
coherence lengths agreed with the HWHM plot better than the least squares 
techniques. The HWHM and iterative fit proved to be the most stable 
parameterizations of the E-field coherence length, and were used in all 
subsequent coherence length comparisons and plots. 

Other factors were also considered in the coherence length calculation. 
As stated earlier, some choices for point source, such as the Airy-type source 
required more energy at high spatial frequencies than others, such as a narrow 
Gaussian. While this difference produced < 2% effect in the normalized 
irradiance variance calculations. It appeared to have more effect on coherence 
length calculations. Figure (53) shows that an Airy-type source produced 
coherence lengths up to 10% higher at low turbulence strengths than those 
from the Gaussian-type source of Eq. (109). Therefore, only the latter type 
source was used for further investigations. 

For an arbitrary spectrum of refractive index fluctuations <I>„(k,z), the 
integral formulation Eq. (61) provided tiie wave structure function for the 
atmospheric MTF in Eq. (134). Spedfically, numerical integration of Eq. (61) for 
the spectra with grid cutoff, Gaussian, and Hill/Frehlich viscous convective 
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Figura 53 Iterative fit coherence lengths versus turbulence strength for Airy- 
type source (dotted), Gaussian-type source (solid), and Gaussian-type source 
with zero Karhunen-Loeve low spatial frequency corrections (dashed). 





enhancement inner scales allowed comparison of theory (Eq. (134)) and the 
simulation (described in Chapter IV). However, in the structure 

function numerical integration, care had to be exercised in properly treating the 
low spatial frequency portion of the integrand since the spectra contained a 
Singularity as k approached zero. This low frequent portion was critical 
to obtain coherence lengths that approached the Kolmogorov theoretical values 
in the Rytov regime. The integral was successfully evaluated by integrating 
analytically for 0 < k < ic^ . (where 1x10*^ m*') (Walters, 1994) since the 
inner scale function F(k^) - 1 here and this portion of the spectrum remains 
Kolmogorov. The remaining portion of the integral that contained the inner 
scale contribution was carried out numericaily. More realistic spectra could also 
have included an outer scale, but again no universal form of outer scale exists 
due to anisotropy of the atmosphere at large scale sizes. When included in the 
numerical integrations for test purposes, an outer scale raised the coherence 
length compared to the Kolmogorov case. However, these investigations did 
not use an explicit outer scale in the simulations. 
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IV. RESULTS 


A. NORMAUZED IRRADIANCE VARIANCE 

The computer simulation guidelines and considerations discussed in 
Chaplir III were implemented to investigate the behavior of the normalized 
irradiance variance and E-field coherence length in the Rytov and saturation 
regimes for the grid cutoff, Gaussian, and Hill/Frehlich viscous-convective 
enhancement inner scales. Specifically, the simulations S 4 >ply to stratosphenc 
propagation with propagation distance L = 200 km, wavelength X - 500 nm. 
strengths of turbulence = [5x10**. 50], and inner scale sizes [0, 15] cm. 

The simulations used a 1024x1024 grid with grid element size given by Eq. 

(92), an Airy-type source modified to produce a final zero turbulence irradiance 
pattern with edges apodized by a Gaussian (Eq. (109)) and width corresponding 
to half the grid width, 32 phase screens utilizing a five-term Karhunen-Loeve 
low spatial frequency correction, and 30 realizations in each set of runs. The 
central 256x256 portion of each propagated E-field was used for the normalized 
irradiance variance and coherence length calculations. 

For the Gaussian inner scale values of 0 (grid cutoff), 5.10, and 15 cm. 
Fig. (54) plots the normalized irradiance variance versus turbulence strength 
in the Rytov regime from both numedcal integration of the equation from Rytov- 
Tatarski theory. Eq. (34) (dotted lines), and from computer simulations 
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Figim M Normalized irradiance variance (divided by from simulation 
(solid) and numerical integration (dotted) for Gaussian inner scales of 0, 5, 10, 
15 cm. 
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(solid tines). All values are normalized by . Numerical integration values 
for 0 cm inner scale agreed within 1% of the theoretical zero inner scale values. 
The difference resulted because the numerical Integration was limited to spatial 
frequencies below the grid cutoff k^. Larger grid cutoff values gave closer 
agreement. The simulation normalized itradiance variances agreed within 2% 
of the numerical integration values for all four inner scale values examined. 
Nonzero Gaussian inner scales reduced ttie normalized irradiance variance 
below the zero inner scale value (by 10%. 25%, and 40% for the 5, 10, 15 cm 
cases, respectively). Intuitively, the finite inner scale suppressed the higher 
spatial frequency index of refraction fluctuations and thus reduced the variance. 

Figure (55) plots the normalized irradiance variance (divided by Pg^) for 
the single turbulence strength - 5x10** and Gaussian Inner scale sizes of 0 
(grid cutoff), 5. 10, and 15 cm. The numerical integration values (dotted line) 
and computer simulation values (solid line) showed an almost linear decrease 
of the normalized irradiance variance with increasing inner scale size in the 
Rytov regime. As Flatty, Wang, and Martin (1993) point out, the Gaussian 
inner scale does not accurately describe ttie inner scale observed in the 
atmosphere but retains usefulness because it facilitates some theoretical 
calculations. 

For the Hill viscous-convective enhancement inner scale sizes of 0 (grid 
cutoff), 5,10, and 15 cm. Fig. (56) plots the normalized irradiance variance in 
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Flgun 65 Normalized irradianoe variance from simulation (solid) and numerical 
integration (dotted) for Gaussian inner scales of 0,5,10.15 cm at s 5x10*^. 
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Figure 66 Normalized irradiance variance (divided by from simulation 
(solid) and numerical integration (dotted) for Hill viscous-convective 
enhancement inner scales of 0,5,10,15 cm. 










r 


the Rytov regime versus turbulence strength . Numerical integration of the 
Rytov-Tatarski theory, Eq. (34) (clotted lines), and computer simulation (solid 
lines) agreed within 2%. For the smaller inner scales, the normalized irradiance 
variance exceeded the zero inner scale values (by 30% for the 5 cm case). 

The viscous-convective enhancement of tt)e strength of higher spatial 
wavenumber fluctuations near the inner scale wavenumber increased the 
variance. Yet, for large enough inner scale, the rolloff beyond the enhancement 
suppressed the higher spatial frequency fluctuations enough to eventually 
reduce the variance below the zero inner scale variance (by 30% for the 15 cm 
case). The 10 cm values happened to lie within 1% of the 0 cm values. 

Figure (57) plots the normalized irradiance variance (divided by for 
the single turbulence strength = 5x10** and the Hill viscous-convective 
enhancement inner scale sizes of 0 (gnd cutoff), 2, 3. 4, 5, 6, 7, 10. and 15 cm. 
Numerical integration of the Rytov-Tatarski results (dotted line) and computer 
simulation (solid line) dearly illustrated the rising and then falling behavior of 
the normalized irradiance variance wifri increasing inner scale size in the Rytov 
regime. The normalized irradiance variance achieved a maximum for ^ ~ 4 cm, 
which was about 30% of the Fresnel length R, ■ l / i n - 12.6 cm. 

The dashed line in Fig. (57) diows simulation values using the Frehlich 
parameterization of the viscous-convective enhancement inner scale. The 
Frehlich inner scale shifted the plot slightly to smaller inner scale sizes, 
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maximized at 2% less than the Hill maximum, and matched the Hill values 
within 3% over the range of inner scale plotted. Additionally, simulation runs 
using 4 cm Hill and Frehlich inner scales agreed within 3% over the range of 
turbulence strengths 5x10** < < 50 . Thus, the Hill and Frehlich versions 

of the viscous-convective enhancement inner scale perform almost identically. 

Previous investigations have illustrated the dramatic monotonic rise in 
normalized irradiance variance in the saturation regime as the inner scale size 
increases (Martin and Flatty, 1988). Figure (58) shows normalized irradiance 
variance from computer simulations for 0 (grid cutoff), 5, 10, and 15 cm 
Gaussian inner scales and a propagation path of 200 km. Figure (59) shows a 
corresponding plot for 0 (grid cutoff). 5.10, and 15 cm Hill inner scales. The 
turbulence values range from the Rytov regime (low turbulence with 1) to 
the saturation regime (high turbulence, sr 1 ). The Rytov regime showed 
again the behaviors illustrated with Figs. (54) - (57). Increasing Gaussian inner 
scale size produced monotonicaliy decreasing normalized irradiance variance. 
Martin and Flatt6 (1988, 1990) provided similar plots of normalized irradiance 
variance with a Gaussian inner scale. Tbe Hill viscous-convective 
enhancement caused the normalized irradiance variance to rise and then fall as 
the inner scale size Increased. However, in the saturation regime, the 
normalized irradiance variance increased monotonicaliy with increasing inner 
scale size for both the Gaussian and Hill inner scales. The transition in 
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Flgura 58 Normalized irradiance variance over Rytov and saturation regimes 
for Gaussian inner scales of 0 (solid), 5 (dashed), 10 (dash-dot), and 15 cm 
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Figure 69 Normalized irradiance variance over Rytov and saturation regimes 
for Hill viscous-convective enhancement inner scales of 0 (solid), 5 (dashed), 
10 (dash-dot), and 15 cm (dotted). 
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behavior occurred with the onset of saturation around ~ 1 • This crossing 
behavior has been plotted for the log intensity variance with the viscous- 
convective inner scale by Hill and Clifford (1978). 

Figures (54) - (57) illustrate the close agreement between numerical 
integration and computer simulation values for the normalized irradiance 
variance at low strengths of turbulence. This agreement provided a validity 
check on these computer simulations that incorporated an inner scale. 

B. COHERENCE LENGTH 

Coherence lengths of the E-field were calculated from the same 
1024x1024 simulation runs used for the normalized irradiance variance 
calculation. The average atmospheric MTF was formed from 30 realizations 
using the FFT autocorrelation method, and then the corresponding coherence 
length rg and HWHM were calculated. The iterative fit rg's were used for 
comparisons because they most closely followed the HWHM behavior. The 
HWHM may provide the coarsest measure of coherence length but, since it 
requires no assumptions about the form of the MTF, it may also be the most 
reliable. 

Figure (60) shows the coherence length Tq from numerical integration of 
Eq. (61) for an approximately zero inner scale and normalized by the theoretical 
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Figure 60 Coherence length from numerical integration for zero inner scale 
versus turbulence strength. Values normalized by Kolmogorov turbulence zero 
inner scale coherence length. 


148 
















Kolmogorov turbulence zero inner scale coherence length rg of Eq. (66). This 
plot indicates that the numerical integration coherence lengms calculated in 
these investigations were accurate within 2%. 

Figure (61) shows the coherence length rg from numerical integration of 
Eq. (61) with the grid cutoff = 318 rad/m for the 1024x1024 grid. The grid 
cutoff did not affect the coherence lengtt) until « 1.5. Above that level of 
turbulence, the absence of the very high spatial frequency contribution to the 
integral caused the coherence length to become larger than the theoretical zero 
inner scale value. The grid cutoff inner scale implicit in the computer 
simulations with the finite grid should have a similar effect at these high 
turbulence values. 

Figure (62) shows the coherence length rg, normalized by the theoretical 
coherence length, from computer simulation of a spherically diverging E-field. 

In the Rytov regime, the simulation agreed with theory to within the ~5% 
overestimation due to using only five terms in the Karhunen-Loeve low spatial 
frequency correction to the phase screens (see Fig. (50)). In the saturation 
rr ime, the simulation coherence length (solid line) dropped below the 
theoretical prediction (by ~25% at Pg^ = 50) and the HWHM (dotted line) mirrored 
this decrease. For strong turbulence, the first order perturbation theory basis 
for coherence length appears to lose validity, as it did for the normalized 
irradiance variance in the saturation regime. 
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FIgura Coherence length from numerical integration for grid cutoff inner 
scale versus turbulence strength. Values normalized by Kolmogorov turbulence 
zero inner scale coherence length. 
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Flgura 62 Coherence length r, (solid line) and HWHM (dotted line) for 
spherical wave propagation through turbulence (Values normalized by 
theoretical coherence lengths). 






The drop in coherence length presumably came from strong scattering in 
high turbulence, but a possible origin in the computer simulation was also 
considered. The simulations started witr. an approximate point source to 
emulate spherical divergence of the E-fieid and required the E-field to remain 
basically confined within the simulation grid over the propagation. The resulting 
E-field properties could have differed from those of a true spherical wave. To 
investigate this possibility, the width of tfie final irradiance pattern was varied 
such that wider final irradiance fields (hence narrower sources) more dosely 
approximated a true point source. Figure (63) plots the results and shows that 
increasing the final irradiance width (diamond, then triangle, then square, then 
X) actually lowered the coherence length in the saturation regime while still 
following the theory in the Rytov regime. This behavior indicated that the 
observed decrease in the saturation regime was physical and not due to 
simulation constraints. 

To further investigate this phenomenon, beam wave (i.e. divergence 
intermediate between spherical and plane waves) and plane wave 
approximations were propagated on a 512x512 grid in which the width of the 
source varied from 4 grid elements (Ax) to 364 grid elements (3/4 of the grid 
width). Figure (64) plots the resulting coherence lengths and indicates that, 
based on the behavior in the Rytov regime, the 4 Ax source (circles) produced 
the spherical wave having large divergence of the E-field, the intermediate 
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FIgm 63 Coherence lengths Fq for varying amounts of spherical divergence, 
hence final irradiance pattern width: diamond 4/8, triangle » 5/8, square = 6/8, 
X-7/6 grid width. 
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Figim 64 Coherence lengths for increasing source width on a 512x512 grid: 
cirdess4Ax (spherical wave); X's^BAx, squaressIGAx. triangles=32Ax, 
diamonds=64Ax, dotss128Ax. asterisks=256^, plusess384Ax (plane wave). 
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source sizes (X^S Ax, square » 16 Ax, triangle = 32 Ax, diamond » 64 Ax) 
produced beam wave divergences, and wide sources (dot - 126 Ax, asterisk = 
256 Ax, plus = 384 Ax) produced plane waves. 

However, as the turbulence strength approached the saturation regime, 
the behaviors changed. The spherical-type propagation coherence length 
dropped ~15% at = 15, the beam wave propagations actually increased in 
coherence length, and the plane wave propagations first decreased ~5% before 
increasing ~15% at = 15. Some of the unevenness in the beam wave 
coherence lengths occurred because smaller regions of the E-field were used to 
calculated the coherence lengths due to the relatively small size and divergence 
of these waves. 

The cause of these behaviors requires hjrther investigation, but a 
hypothesis can be made. As noted earlier, the spherical wave coherence 
length probably decreased below theory for strong turbulence where the Rytov- 
Tatarski first order perturbation theory was no longer valid. Strong scattering 
may have induced spherical divergence of the beam waves, increasing the 
coherence length, and caused the plane wave approximations to diverge like 
beam waves for high turbulence. 

Recapping the above results, investigations of coherence length via 
computer simulation indicated that first order perturbation theory for coherence 
lengths loses validity in the saturation regime, just as it did for normalized 
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irradiance variance. The behavior of coherence length in the saturation regime 
for spherical, beam, and plane wave cases requires further research. 

These investigations then examined the effect of inner scale upon 
coherence length for a spherically diverging beam, utilizing the 1024x1024 
realizations run for the normalized irradiance variance calculations. Figure (65) 
shows the coherence length r^ from numerical integration of Eq. (61) for 
Gaussian inner scales of 5, 10, and 15 cm. The Gaussian inner scale made 
the coherence lengths larger by reducing the energy at high spatial frequencies. 
Larger inner scales monotonically produced larger coherence lengths from 
numerical integration at a given turbulence strength (~8% larger at ~ 1 
(beginning of saturation regime) for ^ = 15 cm). A plot of HWHM from 
numerical integration for Gaussian inner scales would appear similar since the 
theoretical HWHM is proportional to the Kolmogorov coherence length r,, 

HWHM = 0.382 Tq. (147) 

Figure (66) shows the coherence length r^ and Fig. (67) shows the 
HWHM from wave optics computer simulations with Gaussian inner scales of 5, 
10, and 15 cm (solid lines), superimposed upon the numerical integration 
predictions of Fig. (65) (dotted lines). Itie rg and HWHM simulation values 
agreed well with theory for < 0.1, but started deviating frc..i theory even 
before Po^ = 1, i.e. the onset of saturation for the normalized irradiance 
variance. In the saturation regime, computer simulation r^ and HWHM still 
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Figm 66 Coherence length r, from numerical integration for Gaussian inner 
scales of 5, 10, and 15 cm, normalized by the theoretical zero inner scale 
coherence lengths. 
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Figura 66 Coherence lengths Tq from computer simulation with Gaussian inner 
scales of 5, 10, and 15 cm; values normalized by the theoretical zero inner 
scale coherence lengths. 


158 









2.0 


i 1.5 

I 


a 

u 


c 

c 

§ 


S 

X 


0.5 


■ T-i-i m il l I r I Him I I if nwy ^'Tmnn " " i~ri ii ii i j > ii f i n i j 

dotted lines * num. integration HWfHM, GAUSSIAN inner scale 
solid lines » simulotion HWHM. GAUSSIAN inner scole 


1.0^- 





[. pluses 
diomonds 
triongles 
squares 


> inner scole of 0 cm 
«• inner scole of 5 cm 
inner scale of 10 cm 
■ inner scole of 15 cm 


I i tut 


iiL 


0.0 

0.0001 0.0010 


I I III 


luL 


» «»»« 


uiL 


■ I i m 


luL 


■I I > II 


iuL 


I I 11 mil 


0.0100 0.1000 1.0000 10.0000 100.0000 

Beta squored 


Figure 67 HWHM's from computer simulation with Gaussian inner scales of 5, 
10, and 15 cm; values normalized by the theoretical zero inner scale HWHM's. 
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agreed well with each other, and their behavior showed a combination of 
decrease in coherence length due to saturation regime turbulence and inaease 
in coherence length due to inner scale. 

To elucidate the effects of inner scale by itself, Figs. (68) and (69) show 
the same plots of computer simulation coherence length and HWHM, but now 
normalized by the computer simulation zero inner scale and HWHM values, 
respectively, effectively removiny ..j>aturation contribution. These plots 
clearly show that, even in the saturation regime, inner scale increased the 
coherence length (solid lines) similarly to the predictions of theory (dotted lines). 

The next set of figures plots coherence length and HWHM for 
propagations through turbulence with the Hill viscous-convective enhancement 
inner scale. Figure (70) plots the predicted coherence lengths from numerical 
integration of Eq. (61). Note that the numerical integration coherence lengths 
dropped ~5% in the range = [0.1, 1] before increasing for higher turbulence 
strength. Figures (71) and (72) plot the computer simulation coherence lengths 
Tq and HWHM, normalized by the theoretical values for spherical waves. In the 
Rytov regime, the coherence length r^ decreased (~5% for ^ = 15 cm) even for 
very low turbulence (Og^ < 0.1 )z, but then increased in the saturation regime as 
for the Gaussian inner scale. Figures (73) and (74) plot the same coherence 
lengths r^ and HWHM's but divided by the zero inner scale computer simulation 
values to emphasize the effect of inner scale. Figure (73) clearly shows the 
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small decrease in coherence length r, in the Rytov regime, and both plots show 
the general agreement between theory and computer simulation for the effect of 
inner scale upon coherence length in the saturation regime. In summary, the 
inner scale increased the coherence length in the saturation regime as much as 
50% compared to the zero inner scale case, and more than compensated for 
the decrease from strong scatter. 
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Figure 68 Coherence lengths from computer simulation for Gaussian inner 
scales of 5, 10, and 15 cm; values normalized by computer simulation zero 
inner scale rg. 
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Figure 69 HWHM's from computer simulation with Gaussian inner scales of 5, 
10, and 15 cm; valued normalized by computer simulation zero inner scale 
HWHM’s. 
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Figim 70 Coherence length r, from numerical integration for Hill viscous- 
convective enhancement inner scales of 5,10, and 15 cm, normalized by the 
theoretical zero inner scale coherence length. 


164 









FIgtm 71 Coherence lengths r, from computer simulation with Hill viscous- 
convective enhancement inner scales of 5. 10, and 15 cm; values normalized 
by the theoretical zero inner scale coherence lengths. 
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Figura 72 HWHM's from computer simulation with Hill viscous-convective 
enhancement inner scales of 5,10, and 15 cm; values normalized by the 
theoretical zero inner scale HWHM's. 
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Flgiire 73 Coherence lengths from computer simulation for Hill viscous- 
convective enhancement inner scales of 5.10, and 15 cm; values normalized 
by computer simulation zero inner scale rg. 
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Figure 74 HWHM's from computer Emulation with Hill viscous-convective 
enhancement inner scales of 5, 10, and 15 cm; valued normalized by computer 
simulation zero inner scale HWHM's. 
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V. CONCLUSIONS 


Variations in the index of refraction within a turbulent medium alter an 
E-field propagating through the medium. The Rytov-Tatarski, perturbation 
theory predicts the effects of the turbulence upon the irradiance statistics and 
coherence length of the propagating E-field. Computer simulations modeled the 
propagation of the E-field through the turbulent medium, producing irradiance 
and coherence statistics to compare with theoretical results. 

These investigations used a split-step. Huygens-Fresnei, wave optics, 
computer simulation to model an E-field propagating through a turbulent 
stratosphere. The limits of validity of the simulations were determined based 
upon aliasing considerations, choice of source, and robustness of statistical 
calculations, and produced the following guidelines; 

• The element size for an NxN grid should satisfy Ax = l / K . 

• The maximum turbulence strength for which an NxN grid produces valid 
E-fields is given by * 0.1 N . 

• A coherence length rg» 2.5 Ax corresponds to the maximum turbulence 
strength for which an NxN grid produces valid E-fields. 

• The number of phase screens should be ^ 30 . 

• The number of realizations should be ^ 30 . 

• Low spatial frequency corrections to phase screens improve the accuracy 
of the normalized irradiance variance by 5% and of the E-field coherence 
length by 30%. 
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• Half width at half maximum of the atmospheric MTF and an iterative fit 
provide the most stable parameterizations of the E-field coherence length. 

• Telltale signs of aliasing include a fine-grained irradiance pattern, a boxed 
perimeter of the irradiance pattern, and peaking of energy toward the 
center of the computation grid. 

This research investigated the effect of (1) zero inner scale, (2) Gaussian 
inner scale, (3) Hill's and (4) Frehlich's viscous-convective enhancement inner 
scales, and (5) grid cutoff inner scale on the normalized irradiance variance of a 
spherical wave propagating through a turbulent medium. For the Rytov regime, 
the normalized irradiance variances \Mth grid cutoff, Gaussian, and Hill/Frehlich 
viscous-convective enhancement inner scales were compared to the zero inner 
scale case and to the values from numencal integration of the Rytov-Tatarski 
predictions. For low turbulence strengths, the variances obtained from the 
simulations agreed within 2% of the values from the numerical integrations. 

The grid cutoff inner scale, implicit in discrete grid wave optics computer 
simulations, affected the variance negligibly compared to a true zero inner scale 
at low turbulence strengths with the large 1024x1024 grid. Application of a 
Gaussian inner scale reduced the normalized irradiance variance as much as 
40% for small compared to the zero inner scale case. The more realistic 
Hill viscous-convective enhancement inner scale raised the normalized 
irradiance variance by up to 30% for smaller inner scale values, but for larger 
values of inner scale eventually reduced the variance below the zero inner 


170 








scale value as much as 30%. The latter contrasted with the behavior in the 
saturation regime (f)o^ > 1) where larger inner scales continually enhanced the 
normalized irradiance variance. The Frehlic^ parameterization of the viscous- 
convective enhancement gave normalized irradiance values that agreed within 
3% of the Hill inner scale values over the entire range of turbulence strengths 
investigated. 

The coherence of the E-field was studied by computing the average 
atmospheric MTF from the propagated fields. Parameterizing the MTF with a 
coherence length r,, and half width at half max (HWHM) allowed comparison of 
the coherence length of the E-field with the predicted coherence length from the 
Rytov-Tatarski-Fried theory. In the Rytov regime, simulation coherence lengths 
and HWHM's for spherical and plane wave approximations agreed within 5% of 
the theoretical coherence lengths/HWHM’s for zero inner scale. However, in 
the saturation regime, the spherical wave coherence length decreased as much 
as 25% below the theory. Similar decreases resulted for different widths of the 
final E-field. Beam wave approximations gave coherence lengths that 
increased toward the spherical wave values in the saturation regime, while 
plane wave approximations deviated from ~5% below to ~15% above the theory 
in the saturation regime. Increasing inner scale increased the coherence length 
of a spherically diverging E-field by up to 50% relative to the zero inner scale 
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case in the saturation regime. The amount of the increase agreed with 
numerical predictions from the analytic theory. 

Several avenues for further research exist. Larger grid sizes and finer 
mesh could explore the behavior of normalized irradiance variance and 
coherence length at higher turbulence strengths. In particular, the behavior of 
coherence length for spherical, beam, and plane waves in strong turbulence 
conditions need a more precise description. The atmospheric MTF's from 
simulation and theory could be compared directly, rather than through a 
coherenca length parameterization, and the structure functions could also be 
calculated and compared directly to investigate whether the 2/3 power law 
structure function holds in saturation. Simulations could allow C„^(z) to vary 
along the path and could include an outer scale in the spectrum ot refractive 
index fluctuations to study the effects of large scale anisotropy of atmospheric 
turbulence on normalized irradiance variance and coherence length (possibly 
basing the C„^(z) variations upon high frequency radar data that can resolve the 
large scale variations down to about 300 m). Larger ( > 1 Gbyte RAM) and 
faster computing resources would provide the catalyst for all such further 
investigations of the propagation of an E-field through turbulence. 
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APPENDIX 


Sample listing of YAPS event input file: 


i > ^NT S: 


BLPI ^NA JIQN: 


t 7 ;ciebug flag and random number seed 

0.0 0.053033 ;zenlth angle and r_0 at 0.5 microns 

2 ;number of wavelengths 

0.5e-6 0.5e-6 ;list of wavelengths 

'surf 'beam' 0 0 0 2.531058298 0.0 0.0 0.00 


;surf name, vertex location, clear aperture, 
;super-Gaussian exponent, inner scale 

'surf 'atmosr 0 0 6250 8.0 0.0 0.0 0.00 


;surf name, vertex location, dear aperture, 
;super-Gaussian exponent, inner scale 
'prof 1024 1024 0.009886946 'none' 'dummy' 

jsurface size, grid element size, file flags 
'end' ;end of surface summary 

2 1024 1024 ;number of fields and dimensions 

'times' 0.0 0.090002 ;time initialization 

'thread' 0.0 1.0 1 .propagation start 

finit' 1 0.009886946 2 0 0 0 1 0 0 200000 -i-l 


'apsrf 1 1 'beam' 

'aptou' 1 

'chgfcs' 1 0 0 1.0e+30 
'prop' 1 0 0 200000 
•fldcp' 1 2 
'prop' 2 0 0 0 

'openfl' 7data/davis/c118' 11 
'chgfcs' 2 0 0 200000 
'svflddx' 2 11 385 640 385 640 
fldcp' 1 2 


;initialize field 

;apply aperture profile 

jconvert from amplitude/phase to complex 

jchange focus to apply spherical phase 

;back propagate to create source 

;copy field to use as source later 

ipropagate 

;open field output file 

;remove spherical phase from field 

;save field to output file 

;copy source field to working grid 
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(following three steps repeated 32 times to propagate a distance L ) 

'prop' 2 0 0 193750 ;propagate field 

'mkscm' 6250.0 1.0e-18 193750.0 'atmosi' 

icreate phase screen for Az. Cn2, position 
'apsrf 2 1 'atmosi' :apply phase screen to field 

( save the field ) 

'chgfcs' 2 0 0 200000 ;remove spherical phase from field 

'svfiddx' 2 11 385 640 385 640 ;save field to output file 

fidcp' 1 2 '.copy source field to working grid 

( repeat above propagation 30 times } 


'closefi' 11 idose output file 

'end' ;end simulation 
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